################################################################################### # # Draw empirical distribution funcion of a continuous variable # and tit a normal or lognormal model # # Last review: 17/06/2018 # email: ja@fct.unl.pt # ################################################################################### varlistindex <- 5 varunit <- "(ppm)" nclasses <- 15 fitmodel=TRUE modeltype="lognormal" maxY=1.0 discretizacao <- 40 colourbar <- rgb(0,0,1,0.5) colourline <- rgb(1,0,0,0.5) ################################################################################### par(mfrow=c(1,2),mar=c(3.5, 3.5, 2, 1), mgp=c(2.25, 0.8, 0)) # # Duplicate gtd -> gtdnew # gtdnew <- data.frame(gtd[,varlistindex]) names(gtdnew)[1]=names(gtd)[varlistindex] # # Frequencies # tablefreq <- data.frame("", nrow=nclasses+1, ncol=4) minvar <- min(gtdnew[,1]) maxvar <- max(gtdnew[,1]) classint <- c(seq(minvar,maxvar,by=(maxvar-minvar)/nclasses)) # # Calculations & graphics # teste <- hist (gtdnew[,1], breaks=classint, las=1, main="Função densidade de probabilidade", xlim=c(minvar,maxvar), ylim=c(0,maxY), col=colourbar, xlab=paste(names(gtdnew)[1],varunit), ylab="f(x)", freq = FALSE, cex.lab=1.0, cex.axis=1.0, cex.main=0.90, cex.sub=0.5) if (fitmodel) { xfit <- seq(min(gtdnew[,1]),max(gtdnew[,1]),length=discretizacao) if (modeltype == "normal") { yfit <- dnorm(xfit,mean=mean(gtdnew[,1]),sd=sd(gtdnew[,1])) } else if (modeltype == "lognormal") { mmm=log(mean(gtdnew[,1])) ddd=log(1+((sd(gtdnew[,1])^2/(mean(gtdnew[,1])^2)))) mm=mmm-0.5*ddd dd=sqrt(ddd) yfit <- dlnorm(xfit,mean=mm,sd=dd) } lines(xfit, yfit, col=colourline, lwd=2) } a <- c(paste(teste$breaks[1:nclasses],";",teste$breaks[1:nclasses+1],sep="")) b <- teste$counts c <- teste$counts/sum(teste$counts) d <- cumsum(c) tablefreq <- data.frame(a,b,c,d) names(tablefreq) <- c("Intervalo","Freq_absolutas","Freq_relativas","Freq_cumulativas") plot (ecdf(gtdnew[,1]), do.points=FALSE, verticals=TRUE, main="Função de distribuição cumulativa", xlab=paste(names(gtdnew)[1],varunit), ylab="F(x)", cex.lab=1.0, cex.axis=1.0, cex.main=0.90, cex.sub=0.5) if (fitmodel) { xfit <- seq(min(gtdnew[,1]),max(gtdnew[,1]),length=discretizacao) if (modeltype == "normal") { yfit <- pnorm(xfit,mean=mean(gtdnew[,1]),sd=sd(gtdnew[,1])) text (minvar+0.5*(maxvar-minvar),0.05,paste("Normal (m=",round(mean(gtdnew[,1]),digits=3),"sd=",round(sd(gtdnew[,1]),digits=2),")")) } else if (modeltype == "lognormal") { yfit <- plnorm(xfit,mean=mm,sd=dd) text (minvar+0.5*(maxvar-minvar),0.05,paste("Logormal (m=",round(mm,digits=3),"sd=",round(dd,digits=2),")")) } lines(xfit, yfit, col=colourline, lwd=2) }