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)