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)