A explicação do script a seguir está no artigo R na Prática – Parte 25.
hquebras=function(x){
n=length(x)
nc=ceiling(1+3.222*log10(n))
zc=(max(x)-min(x))/nc
quebras=c(rep(0,(nc+1)))
for (i in 1:(nc+1)){
quebras[i]=min(x)+(i-1)*zc}
return(quebras)
}
setwd("C:\\geoStats\\dados\\jura_goovaerts")
dados=read.csv("jura_prediction.CSV",sep=";",
header=TRUE)
#setwd("C:\\geoStats\\figuras\\chap3")
#jpeg("Fig3_13.jpeg",width=5,height=5,
#units="in",res=300)
par_atuais<-par(no.readonly=TRUE)
par(mfrow=c(2,2))
hist(dados$Cd,col="yellow",xlab="Cd (ppm)",
ylab="Frequência",breaks=hquebras(dados$Cd),
main="Distribuição de Cd (ppm)")
hist(dados$Pb,col="red",breaks=hquebras(dados$Pb),
main="Distribuição de Pb (ppm)",
xlab="Pb (ppm)",ylab="Frequência")
hist(dados$Zn,col="blue",
main="Distribuição de Zn (ppm)",
xlab="Zn (ppm)",ylab="Frequência")
hist(dados$Ni,col="green",
main="Distribuição de Ni (ppm)",
xlab="Ni (ppm)",ylab="Frequência")
par(par_atuais)
#dev.off()