source('TailInflationB.R') # set.seed(14) # ----------------------------------------------------------- # # Basic example for a piecewise linear theta: # Generate data from P such that # P(dx) = e^theta(x) P0(dx): n <- 400 alpha <- 1 beta <- 5 # Preliminary choice of tau and theta for rate beta = 1: t.true <- c(0, 2.0, 4.0, 6.0) th.true <- c(0, 0.5, 1.5, 2.7, 0.8) # Normalization of theta: th.true <- LocalNormalize.2B(t.true, th.true, alpha) # Rescaling to actual beta: t.true <- t.true/beta th.true <- c(th.true[1:length(t.true)], th.true[1+length(t.true)]*beta) # Generate data x <- Simulate.2B(n, t.true, th.true, alpha, beta) 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.B0(x,w,alpha,beta) # Same computation without graphics: res <- TailInflation.B(x,w,alpha,beta) t.hat <- res\$tau th.hat <- res\$theta tt <- seq(0,max(x),length.out=501) tt <- unique(sort(c(tt,t.hat))) # Depicting directional derivatives to check # optimality: x0 <- x*beta t.hat0 <- t.hat*beta th.hat0 <- c(th.hat[1:length(t.hat)], th.hat[1+length(t.hat)]/beta) tt0 <- tt*beta H.tt0 <- LocalDirDeriv2.2B(tt0,x0,w,t.hat0,th.hat0,alpha) plot(tt0,H.tt0[,1],type='l', xlab=expression(italic(t)), ylab=expression(italic(h(t)))) abline(h=0,col='red') abline(v=t.hat0,col='blue') rug(x0) # The function h should be <= 0 with equality at t.hat... # Display true (green) and estimated (black) # tail inflation function: th.true.tt <- Evaluate.2B(tt,t.true,th.true) D.tt <- LocalLinearSplines1.2B(t.hat,tt) th.hat.tt <- D.tt %*% th.hat plot(tt,th.true.tt,type='l',lwd=2,col='green', xlab=expression(italic(x)), ylab=expression(italic(theta(x))), ylim=range(c(th.true.tt,th.hat.tt))) lines(tt,th.hat.tt,lwd=2) rug(x) # Display true (green), estimated (black) and # reference density (magenta): p0.tt <- dgamma(tt,alpha,beta) p.true.tt <- exp(th.true.tt)*p0.tt p.hat.tt <- exp(th.hat.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(p0.tt[2],p.true.tt[2],p.hat.tt[2]))) lines(tt,p.true.tt,col='green',lwd=2) lines(tt,p.hat.tt,lwd=2) rug(x) # ----------------------------------------------------------- # # Another example: scale mixtures of chi_1^2 r.v.s: n1 <- 200 n2 <- 150 n3 <- 50 S2 <- 1.5 S3 <- 2 n <- n1+n2+n3 w <- rep(1/n,n) x <- rchisq(n,df=1)*c(rep(1,n1),rep(S2,n2),rep(S3,n3)) alpha <- 1/2 beta <- 1/2 # 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.B0(x,w,alpha,beta) # Same computation without graphics: res <- TailInflation.B(x,w,alpha,beta) t.hat <- res\$tau th.hat <- res\$theta tt <- seq(0,max(x),length.out=501) tt <- unique(sort(c(tt,t.hat))) # Depicting directional derivatives to check # optimality: x0 <- x*beta t.hat0 <- t.hat*beta th.hat0 <- c(th.hat[1:length(t.hat)], th.hat[1+length(t.hat)]/beta) tt0 <- tt*beta H.tt0 <- LocalDirDeriv2.2B(tt0,x0,w,t.hat0,th.hat0,alpha) plot(tt0,H.tt0[,1],type='l', xlab=expression(italic(t)), ylab=expression(italic(h(t)))) abline(h=0,col='red') abline(v=t.hat0,col='blue') rug(x0) # The function h should be <= 0 with equality at t.hat... # Display true (green) and estimated (black) # tail inflation function: th.true.tt <- log((n1 + n2*S2^(-alpha)*exp((beta - beta/S2)*tt) + n3*S3^(-alpha)*exp((beta - beta/S3)*tt))/n) D.tt <- LocalLinearSplines1.2B(t.hat,tt) th.hat.tt <- D.tt %*% th.hat plot(tt,th.true.tt,type='l',lwd=2,col='green', xlab=expression(italic(x)), ylab=expression(italic(theta(x))), ylim=range(c(th.true.tt,th.hat.tt))) lines(tt,th.hat.tt,lwd=2) rug(x) # Display true (green), estimated (black) and # reference density (magenta): p0.tt <- dgamma(tt,alpha,rate=beta) p.true.tt <- (n1*dgamma(tt,alpha,rate=beta) + n2*dgamma(tt,alpha,rate=beta/S2) + n3*dgamma(tt,alpha,rate=beta/S3))/n p.hat.tt <- exp(th.hat.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(p0.tt[2],p.true.tt[2],p.hat.tt[2]))) lines(tt,p.true.tt,col='green',lwd=2) lines(tt,p.hat.tt,lwd=2) rug(x) # Zoom in by a factor of gamma in x-direction: gamma <- 10 plot(tt,p0.tt,type='l',lwd=2,col='magenta', xlab=expression(italic(x)), ylab=expression(italic(p(x))), xlim=c(0,max(tt)/gamma), ylim=c(0,max(p.true.tt[2],p.hat.tt[2]))) lines(tt,p.true.tt,col='green',lwd=2) lines(tt,p.hat.tt,lwd=2) rug(x)