################################################################################### # # Generate and view a lognormal distribution # # Last review: 17/06/2018 # email: ja@fct.unl.pt # ################################################################################### ################################################################################### opcao=1 # or 2 media=2 desvp=2 ndecimals <- 2 nvalores <- 1000000 discretizacao <- 200 colour <- 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)) # # Calculations # if (opcao == 1) { variancia=desvp^2 medlog=exp(media+0.5*desvp^2) varlog=exp(2*media+desvp^2)*(exp(desvp^2)-1) } else { medlog=media varlog=desvp^2 media=log(medlog)-0.5*log(1+(varlog/medlog^2)) variancia=log((varlog/medlog^2)+1) desvp=sqrt(variancia) } xmin=0 xmax=2*medlog xmaxmin=xmax-xmin tolx=0.1*xmaxmin x=seq(xmin,xmax,length=discretizacao) ydist=dlnorm(x,mean=media,sd=desvp) yprob=plnorm(x,mean=media,sd=desvp) valores=rlnorm(nvalores,mean=media,sd=desvp) # # Graphics # media=round(10^ndecimals*media)/(10^ndecimals) desvp=round(10^ndecimals*desvp)/(10^ndecimals) medlog=round(10^ndecimals*medlog)/(10^ndecimals) varlog=round(10^ndecimals*varlog)/(10^ndecimals) plot(x,ydist,type="l",xlim=c(min(x)-tolx,max(x)+tolx),ylim=c(0,1.5*max(ydist)), lwd=2,col="red",xlab="x",ylab="f(x)", main = paste("Função probabilidade\nLei lognormal média=",media,", desvio padrão=",desvp), cex.lab=1.0, cex.axis=1.0, cex.main=0.90, cex.sub=0.5) legend(xmin, 1.5*max(ydist), legend=c(paste("E(X)=",round(medlog,digits=ndecimals)), paste("VAR(X)=",round(varlog,digits=ndecimals)), paste("DESVP(X)=",round(sqrt(varlog),digits=ndecimals))),cex=0.75) polygon(c(xmin,x,xmax,xmin),c(0,ydist,0,0),col="lightgray",border=NA) lines(x,ydist,type="l",lwd=2,col=colour) points(medlog,0,pch=8,col=rgb(0,0,1,0.5)) plot(x,yprob,type="l",xlim=c(min(x)-tolx,max(x)+tolx),ylim=c(0,1), lwd=2,col="red",xlab="x",ylab="F(x)",main="Função distribuição cumulativa", cex.lab=1.0, cex.axis=1.0, cex.main=0.90, cex.sub=0.5) polygon(c(xmin,x,xmax,xmin),c(0,yprob,0,0),col="lightgray",border=NA) lines(x,yprob,type="l",lwd=2,col=colour) # # Results # print ("Lognormal distribution law simulated values: ") print.table(round(valores, digits=ndecimals)) mediana=qlnorm(0.5, media, desvp) quantil=qlnorm(pnorm(media+desvp, mean=media, sd=desvp), media, desvp) log(mediana) log(quantil)-media print (paste("Média:", round(media, digits=ndecimals), round(medlog, digits=ndecimals))) print (paste("Variância:", round(variancia, digits=ndecimals), round(varlog, digits=ndecimals)))