source('LogConDens.R') # Simulate data: # Sample size: n <- 500 # # Standard Gaussian distribution: # x <- sort(rnorm(n)) # # Standard exponential distribution: # x <- sort(-log(runif(n))) # # Standard logistic distribution: # u <- sort(runif(n)) # x <- log(u/(1-u)) # # Gamma distribution: # shape <- 3 # x <- sort(rgamma(n,shape)) # Beta distribution: shape1 <- 2 shape2 <- 3 x <- sort(rbeta(n,shape1,shape2)) # Estimation of concave log-density, first with # graphical displays of all intermediate steps # (make sure that the graphics window is visible beside the # console window, and the latter should be in the foreground): res <- LogConDens0(x,m0=1) res <- LogConDens0(x,m0=6) m1 <- floor(sqrt(n)) res <- LogConDens0(x,m0=m1) # m0 + 1 is the number of knots for the very first local # search, starting from a Gaussian fit. # The same algorithm without graphical displays: system.time(res <- LogConDens(x,m0=1)) system.time(res <- LogConDens(x,m0=2)) system.time(res <- LogConDens(x,m0=6)) system.time(res <- LogConDens(x,m0=m1)) system.time(res <- LogConDens(x,m0=2*m1)) system.time(res <- LogConDens(x,m0=3*m1)) # Very often, the choice m0 = sqrt(n) yields the # shortest running time. # Plot of directional derivatives for additional # knots (should all be <= 0 with equality at the # actual estimated knots): par(mai=c(0.8,0.9,0.1,0.1)) plot(x,res\$DLtau,type='l', xlab=expression(italic(tau)), ylab=expression(italic(DL(phi,V[tau])))) rug(x) abline(h=0,col='red') points(res\$knots,rep(0,length(res\$knots)),col='blue') # Plot of estimated (black) and true (green) log-density: xx <- LocalInterpolate.1A(x,7) phixx <- LocalInterpolate.1A(res\$phix,7) # # Standard Gaussian truth: # f0xx <- dnorm(xx) # # Standard logistic truth: # f0xx <- 1/(exp(xx) + exp(-xx) + 2) # phi0xx <- log(f0xx) # # Gamma truth: # f0xx <- dgamma(xx,shape) # phi0xx <- log(f0xx) # Beta truth: f0xx <- dbeta(xx,shape1,shape2) phi0xx <- log(f0xx) plot(xx,phi0xx,type='l',lwd=3,col='green', xlab=expression(italic(x)), ylab=expression(italic(phi(x))), ylim=range(c(res\$phi,phi0xx))) lines(xx,phi0xx,col='darkgreen') lines(res\$x,res\$phix,lwd=2) points(res\$knots,res\$phi,pch=16) rug(x) # Plot of estimated (black) and true (green) density: plot(xx,f0xx,type='l',lwd=3,col='green', xlab=expression(italic(x)), ylab=expression(italic(f(x))), ylim=c(0,max(c(exp(res\$phi),f0xx)))) lines(xx,f0xx,col='darkgreen') lines(xx,exp(phixx),lwd=2) rug(x)