source('TailInflationA.R')

# Simulate data from N(mu,sigma^2):
mu <- 0.5
sigma <- 1.25
n <- 400
x <- sort(rnorm(n,mu,sigma))
w <- rep(1/n,n)

# Estimation of tail inflation function with
# graphical display of all intermediate steps
# (make sure the console window is in the foreground
#  and graphics window beside but visible):
res <- TailInflation.A0(x,w)
# Same computation without graphics:
res <- TailInflation.A(x,w)
res
tau <- res$tau
theta <- res$theta

tt <- seq(min(c(mu-4*sigma,x)),
	max(c(mu+4*sigma,x)),length.out=401)
tt <- unique(sort(c(x,tau,tt)))

# Check optimality of fit by plotting directional
# derivatives for extra kinks:
par(cex=1.5,mai=c(1.0,1.1,0.1,0.05),mgp=c(2.3,1,0))
H.tt <- LocalDirDeriv2.2A(tt,x,w,tau,theta)
plot(tt,H.tt[,1],type='l',
	xlab=expression(italic(t)),
	ylab=expression(italic(h(t))))
abline(h=0,col='red')
for (j in 1:length(tau)){
	abline(v=tau[j],col='blue')
}
rug(x)

# Display true (green) and estimated (black)
# tail inflation function:
th0.tt <- log(dnorm(tt,mu,sigma)/dnorm(tt))
D.tt <- LocalLinearSplines1.2A(tau,tt)
thh.tt <- D.tt %*% theta
plot(tt,th0.tt,type='l',lwd=2,col='green',
	xlab=expression(italic(x)),
	ylab=expression(italic(theta(x))),
	ylim=range(c(th0.tt,thh.tt)))
lines(tt,thh.tt,lwd=2)
rug(x)

# Display true (green), estimated (black) and
# reference density (magenta):
p0.tt <- dnorm(tt)
f0.tt <- exp(th0.tt)*p0.tt
fh.tt <- exp(thh.tt)*p0.tt
plot(tt,p0.tt,type='l',lwd=2,col='magenta',
	xlab=expression(italic(x)),
	ylab=expression(italic(p(x))),
	ylim=c(0,max(c(p0.tt,f0.tt,fh.tt))))
lines(tt,f0.tt,col='green',lwd=2)
lines(tt,fh.tt,lwd=2)
rug(x)


# The same procedures applied to symmetrized data:

x.symm <- sort(c(x,-x))
w.symm <- rep(1/(2*n),2*n)
res <- TailInflation.A0(x.symm,w.symm)
res <- TailInflation.A(x.symm,w.symm)
res
tau <- res$tau
theta <- res$theta

tt <- seq(min(c(-abs(mu)-4*sigma,x,-x)),
	max(c(abs(mu)+4*sigma,x,-x)),length.out=401)
tt <- unique(sort(c(x,-x,tau,tt)))

# Check optimality of fit by plotting directional
# derivatives for extra kinks:
H.tt <- LocalDirDeriv2.2A(tt,x.symm,w.symm,tau,theta)
plot(tt,H.tt[,1],type='l',xlab='t',ylab='h(t)')
abline(h=0,col='red')
for (j in 1:length(tau)){
	abline(v=tau[j],col='blue')
}
rug(x.symm)

# Display true (green) and estimated (black)
# tail inflation function:
th0.tt <- log((dnorm(tt,mu,sigma) + dnorm(tt,-mu,sigma))/
	(2*dnorm(tt)))
D.tt <- LocalLinearSplines1.2A(tau,tt)
thh.tt <- D.tt %*% theta
plot(tt,th0.tt,type='l',lwd=2,col='green',
	xlab=expression(italic(x)),
	ylab=expression(italic(theta(x))),
	ylim=range(c(th0.tt,thh.tt)))
lines(tt,thh.tt,lwd=2)
rug(x.symm)

# Display true (green), estimated (black) and
# reference density (magenta):
p0.tt <- dnorm(tt)
f0.tt <- exp(th0.tt)*p0.tt
# Estimated density:
fh.tt <- exp(thh.tt)*p0.tt
plot(tt,p0.tt,type='l',lwd=2,col='magenta',
	xlab=expression(italic(x)),
	ylab=expression(italic(f(x))),
	ylim=c(0,max(c(p0.tt,f0.tt,fh.tt))))
lines(tt,f0.tt,col='green',lwd=2)
lines(tt,fh.tt,lwd=2)
rug(x.symm)


# Simulate data logistic distribution with
# scale parameter sigma:
sigma <- sqrt(0.5)
n <- 400
u <- runif(n)
x <- sigma*log(u/(1-u))
w <- rep(1/n,n)

# Estimation of tail inflation function with
# graphical display of all intermediate steps
# (make sure the console window is in the foreground
#  and graphics window beside but visible):
res <- TailInflation.A0(x,w)
# Same computation without graphics:
res <- TailInflation.A(x,w)
res
tau <- res$tau
theta <- res$theta

xmax <- max(qnorm(0.9999),log(0.9999/0.0001)*sigma,max(abs(x)))
tt <- seq(-xmax,xmax,length.out=401)
tt <- unique(sort(c(x,tau,tt)))

# Check optimality of fit by plotting directional
# derivatives for extra kinks:
par(cex=1.5,mai=c(1.0,1.1,0.1,0.05),mgp=c(2.3,1,0))
H.tt <- LocalDirDeriv2.2A(tt,x,w,tau,theta)
plot(tt,H.tt[,1],type='l',
	xlab=expression(italic(t)),
	ylab=expression(italic(h(t))))
abline(h=0,col='red')
for (j in 1:length(tau)){
	abline(v=tau[j],col='blue')
}
rug(x)

# Display true (green) and estimated (black)
# tail inflation function:
p.tt <- sigma^(-1)*(exp(tt/sigma) + exp(-tt/sigma) + 2)^(-1)
p0.tt <- dnorm(tt)
th0.tt <- log(p.tt/p0.tt)
D.tt <- LocalLinearSplines1.2A(tau,tt)
thh.tt <- D.tt %*% theta
plot(tt,th0.tt,type='l',lwd=2,col='green',
	xlab=expression(italic(x)),
	ylab=expression(italic(theta(x))),
	ylim=range(c(th0.tt,thh.tt)))
lines(tt,thh.tt,lwd=2)
rug(x)

# Display true (green), estimated (black) and
# reference density (magenta):
ph.tt <- exp(thh.tt)*p0.tt
plot(tt,p0.tt,type='l',lwd=2,col='magenta',
	xlab=expression(italic(x)),
	ylab=expression(italic(p(x))),
	ylim=c(0,max(c(p0.tt,p.tt,ph.tt))))
lines(tt,p.tt,col='green',lwd=2)
lines(tt,ph.tt,lwd=2)
rug(x)