# inst/exploreLogitNormal.R In logitnorm: Functions for the Logitnormal Distribution

```# TODO: Add comment
#
# Author: twutz
###############################################################################

library(ggplot2)

xGrid = c(0+.Machine\$double.eps,seq(0,1, length.out=81)[-c(1,81)],1-.Machine\$double.eps)

#explore chanign sigma at mu=0
theta0 <- expand.grid(mu=seq(0,2,length.out=9), sigma=10^seq(-0.5,0.5,length.out=5))
n <- nrow(theta0)

.calcDensityGrid <- function(
### Calculate lognormal desnity for given combinations
theta0	##<< matrix with columns mu and sigma
,xGrid = seq(0,1, length.out=81)[-c(1,81)]
){
dx <- apply( theta0, 1, function(theta0i){
dx <- dlogitnorm(xGrid, mu=theta0i[1], sigma=theta0i[2])
})
dimnames(dx) <- list(iX=NULL,iTheta=NULL)
ds <- melt(dx)
ds[1:10,]
ds\$mu <- rep(as.factor(round(theta0[,1],2)), each=length(xGrid))
ds\$sigma <- rep(as.factor(round(theta0[,2],2)), each=length(xGrid))
ds\$x <- rep(xGrid, nrow(theta0))
ds
### data frame with columns value,x,mu and sigma
}

.calcCdfGrid <- function(
### Calculate lognormal desnity for given combinations
theta0	##<< matrix with columns mu and sigma
,xGrid = seq(0,1, length.out=81)[-c(1,81)]
){
dx <- apply( theta0, 1, function(theta0i){
dx <- plogitnorm(xGrid, mu=theta0i[1], sigma=theta0i[2])
})
dimnames(dx) <- list(iX=NULL,iTheta=NULL)
ds <- melt(dx)
ds[1:10,]
ds\$mu <- rep(as.factor(round(theta0[,1],2)), each=length(xGrid))
ds\$sigma <- rep(as.factor(round(theta0[,2],2)), each=length(xGrid))
ds\$x <- rep(xGrid, nrow(theta0))
ds
### data frame with columns value,x,mu and sigma
}

#qplot(xGrid,value,data=ds[ds\$mu==0,], geom="line", color=sigma )
windows(width=7,height=4.5)

ds <- .calcDensityGrid(theta0,xGrid=xGrid)
qplot(xGrid,value,data=ds, geom="line", color=sigma, ylab="density")+ facet_wrap(~mu,scales="free")+opts(axis.title.x = theme_blank())
#savePlot("logitnormDensityGrid",type="pdf")

windows(width=4,height=4)
ds <- .calcDensityGrid(theta0,xGrid=xGrid)
qplot(xGrid,value,data=ds[ds\$mu %in% c(0,1),], geom="line", color=sigma, ylab="density")+ facet_grid(mu~.,scales="free")+opts(axis.title.x = theme_blank())
#savePlot("logitnormDensityGrid2",type="pdf")

ds <- .calcCdfGrid(theta0,xGrid=xGrid)
qplot(xGrid,value,data=ds[ds\$mu %in% c(0,1),], geom="line", color=sigma, ylab="cumulative density")+ facet_grid(mu~.,scales="free")+opts(axis.title.x = theme_blank())
#savePlot("logitnormCdfGrid2",type="pdf")

mle=0.8
mu <- seq(0,logit(mle),length.out=10)[-10]
tmp <- twSigmaLogitnorm(mle,mu)
sigma2 <- (logit(mle)-mu)/(2*mle-1)
theta0 <- cbind(mu,sigma=as.numeric(sqrt(sigma2)))

ds <- .calcDensityGrid(theta0,xGrid=xGrid)
qplot(xGrid,value,data=ds, geom="line", color=mu )
qlogitnorm(p=0.975,mu=theta0[,1],sigma=theta0[,2])
qlogitnorm(p=0.99,mu=theta0[,1],sigma=theta0[,2])

mle=0.8
#mtrace(.ofLogitnormMLE)
plot( .ofLogitnormMLE(mu,mle=mle,quant=0.95,perc=0.999) ~ mu)
#mtrace(twCoefLogitnormMLE)
theta <- twCoefLogitnormMLE( mle=mle,quant=0.95,perc=0.99)
plot( dlogitnorm(xGrid,mu=theta[1],sigma=theta[2])~xGrid, type="l")
c( (logit(mle)-theta[1])/(2*mle-1), theta[2]^2 )
abline(v=c(mle,0.99),col="gray")

mle=0.5
#mtrace(twCoefLogitnormMLE)
theta <- twCoefLogitnormMLE( mle=mle,quant=0.99,perc=0.999)
q2 <- qlogitnorm(0.999,mu=theta[1],sigma=theta[2])
plot( dlogitnorm(xGrid,mu=theta[1],sigma=theta[2])~xGrid, type="l")
abline(v=c(mle,q2),col="gray")

mle=0.1
#mtrace(twCoefLogitnormMLE)
theta <- twCoefLogitnormMLE( mle=mle,quant=0.5,perc=0.99)
quant2 <- qlogitnorm(0.99,mu=theta[1],sigma=theta[2])
plot( dlogitnorm(xGrid,mu=theta[1],sigma=theta[2])~xGrid, type="l")
abline(v=c(mle,quant2),col="gray")

mle=0.9
#mtrace(twCoefLogitnormMLE)
theta <- twCoefLogitnormMLE( mle=mle,quant=0.98,perc=0.99)
#theta <- twCoefLogitnormMLE( mle=mle,quant=0.99,perc=0.999)
quant2 <- qlogitnorm(0.99,mu=theta[1],sigma=theta[2])
plot( dlogitnorm(xGrid,mu=theta[1],sigma=theta[2])~xGrid, type="l")
abline(v=c(mle,quant2),col="gray")

#------- sigma for the flattest case
mle=0.9
mu=seq(0,if(mle<0.5) mle else 1-mle,length.out=10)[-10]
mu=seq(0,0.25,length.out=10)[-10]
tmp <- twSigmaLogitnorm(mle,mu)
sigma2 <- (logit(mle)-mu)/(2*mle-1)
theta0 <- cbind(mu,sigma=as.numeric(sqrt(sigma2)))
ds <- .calcDensityGrid(theta0,xGrid=xGrid)
qplot(xGrid,value,data=ds, geom="line", color=mu, ylim=c(0.9,1.2) )
```

## Try the logitnorm package in your browser

Any scripts or data that you put into this service are public.

logitnorm documentation built on May 31, 2017, 4:45 a.m.