### Main programs: TailInflation.B0 <- function(X,W=rep(1/length(X),length(X)), alpha=1,beta=1, xmin = 0, xmax = 1.1*max(X), delta1=10^(-10)/length(X), delta2=10^(-3)*sqrt(alpha)/length(X)) # Computation of a log-convex and increasing density ratio with respect to # the Gamma distribution with shape alpha and rate beta.. # Graphical display of all intermediate steps... { # If X ~ P, where P(dx) = e^theta(x) P_alpha,beta(dx) # with P_alpha,beta the Gamma distribution with # shape alpha and rate beta, then beta*X ~ Pt, where # Pt(dx) = e^thetat(x) P_alpha,1(dx) and thetat(x) = theta(x/beta). # Hence we multiply the data by beta, estimate thetat and obtain # theta(x) = thetat(x*beta) in the end. X <- X*beta # Preparations of raw data: x <- unique(sort(X)) n <- length(x) W <- W/sum(W) w <- rep(0,n) for (i in 1:n) { tmp <- (X == x[i]) w[i] <- sum(W[tmp]) } # MLE for theta in V_{0} \cap \Theta: mu <- sum(w*x) kappa <- max(1 - alpha/mu,0) # \hat{\kappa}^+ ckappa <- -alpha*log(1 - kappa) # c(\hat{\kappa}^+) LL <- -ckappa + kappa*mu - exp(-ckappa)/(1-kappa)^alpha + 1 tmp <- Optimality1.2B(x,w,mu,alpha) if (length(tmp$tnew) > 0){ ht0 <- max(tmp$htnew) t0 <- tmp$tnew[which.max(tmp$htnew)] }else{ ht0 <- 0 } # Graphical display: message <- readline("First plot:") plot(c(xmin,xmax),rep(0,2),type='l', lwd=3,col='green', xlab=expression(italic(x)), ylab=expression(italic(theta(x))), xlim=range(x), ylim=range(kappa*c(xmin,xmax) - ckappa), main=paste('Starting point: ', 'LL = ',round(LL,digits=9),sep='')) lines(c(xmin,xmax),kappa*c(xmin,xmax)-ckappa,lwd=2) rug(x) # Bookkeeping: nr.local.search <- 0 nr.Newton <- 0 if (ht0 < delta2){ print(paste('Parametric fit sufficient: Gamma(', alpha,',',beta,'). Log-likelihood = 0.',sep='')) return(list(tau=NULL,theta=NULL)) } tau <- 0 theta <- c(-ckappa,kappa) # Graphical display: abline(v=tau,col='red') # Check for global optimality: tmp <- Optimality2.2B(x,w,tau,theta,alpha) if (length(tmp$tnew) > 0){ ht0 <- max(tmp$htnew) t0 <- tmp$tnew[which.max(tmp$htnew)] }else{ ht0 <- 0 } while (ht0 > delta2){ nr.local.search <- nr.local.search + 1 # Store current value of LL: LL.old <- LL # Add a new knot t0 to tau: j <- sum(tau <= t0) m <- length(tau) if (j == 0){ # Then tau[1] > 0 tau <- c(t0,tau) theta <- c(theta[1],theta) } if (j > 0 && j < m){ theta1 <- ((tau[j+1] - t0)*theta[j] + (t0 - tau[j])*theta[j+1])/(tau[j+1] - tau[j]) tau <- c(tau[1:j],t0,tau[(j+1):m]) theta <- c(theta[1:j],theta1,theta[(j+1):(m+1)]) } if (j == m){ theta1 <- theta[m] + (t0 - tau[m])*theta[m+1] tau <- c(tau,t0) theta <- c(theta[1:m],theta1,theta[m+1]) } # Graphical display: message <- readline("New knot:") xtx <- c(xmin,tau,xmax) Dtmp <- LocalLinearSplines1.2B(tau,xtx) thetaxtx <- Dtmp %*% theta plot(c(xmin,xmax),rep(0,2),type='l', lwd=3,col='green', xlab=expression(italic(x)), ylab=expression(italic(theta(x))), xlim=range(x), ylim=range(thetaxtx), main=paste('LL =', round(LL,digits=9))) lines(xtx,thetaxtx,lwd=2) points(tau,thetaxtx[2:(length(tau)+1)]) abline(v=tau,col='red') abline(v=t0,col='blue') rug(x) wt <- LocalLinearSplines2.2B(tau,x,w) proposal <- LocalNewton.2B(wt,tau,theta,alpha,delta1) nr.Newton <- nr.Newton + 1 while (proposal$dirderiv > delta1){ LL.Newton <- LL theta.new <- proposal$theta.new # Graphical display: message <- readline("New proposal:") thetaxtx.new <- Dtmp %*% theta.new plot(c(xmin,xmax),rep(0,2),type='l', lwd=3,col='green', xlab=expression(italic(x)), ylab=expression(italic(theta(x))), xlim=range(x), ylim=range(thetaxtx), main=paste('LL =', round(LL.Newton,digits=9))) lines(xtx,thetaxtx,lwd=2) lines(xtx,thetaxtx.new,lwd=2,col='blue') points(tau,thetaxtx.new[2:(length(tau)+1)],col='blue') abline(v=tau,col='red') abline(v=t0,col='blue') rug(x) # 2nd step size correction: corr <- LocalStepSize2.2B(tau,theta,theta.new) if (length(corr$tau) < length(tau)){ tau <- corr$tau theta <- corr$theta wt <- LocalLinearSplines2.2B(tau,x,w) # Graphical display: message <- readline("2nd step size correction:") xtx <- c(xmin,tau,xmax) Dtmp <- LocalLinearSplines1.2B(tau,xtx) thetaxtx <- Dtmp %*% theta lines(xtx,thetaxtx,lwd=2,col='magenta') }else{ theta <- theta.new } # Normalization: theta <- LocalNormalize.2B(tau,theta,alpha) LL <- sum(wt*theta) if (LL > LL.Newton){ proposal <- LocalNewton.2B(wt,tau,theta,alpha,delta1) nr.Newton <- nr.Newton + 1 # For graphical display: thetaxtx <- Dtmp %*% theta }else{ proposal$dirderiv <- 0 } } # Graphical display: message <- readline("Local optimum:") thetaxtx <- Dtmp %*% theta plot(c(xmin,xmax),rep(0,2),type='l', lwd=3,col='green', xlab=expression(italic(x)), ylab=expression(italic(theta(x))), xlim=range(x), ylim=range(thetaxtx), main=paste('LL =', round(LL,digits=9))) lines(xtx,thetaxtx,lwd=2) points(tau,thetaxtx[2:(length(tau)+1)]) abline(v=tau,col='red') rug(x) # Check for global optimality: if (LL > LL.old){ tmp <- Optimality2.2B(x,w,tau,theta,alpha) if (length(tmp$tnew) > 0){ ht0 <- max(tmp$htnew) t0 <- tmp$tnew[which.max(tmp$htnew)] }else{ ht0 <- 0 } }else{ ht0 <- 0 } } # Graphical display: message <- readline("Final estimate:") thetaxtx <- Dtmp %*% theta plot(c(xmin,xmax),rep(0,2),type='l', lwd=3,col='green', xlab=expression(italic(x)), ylab=expression(italic(theta(x))), xlim=range(x), ylim=range(thetaxtx), main=paste('LL =', round(LL,digits=9))) lines(xtx,thetaxtx,lwd=2) points(tau,thetaxtx[2:(length(tau)+1)]) abline(v=tau,col='red') rug(x) # Rescale to account for beta != 1 x <- x/beta tau <- tau/beta theta[length(tau)+1] <- beta*theta[length(tau)+1] return(list(tau=tau,theta=theta, x=x,w=w,LL=LL, nr.local.search=nr.local.search,nr.Newton=nr.Newton)) } TailInflation.B <- function(X,W=rep(1/length(X),length(X)), alpha=1,beta=1, delta1=10^(-10)/length(X), delta2=10^(-3)*sqrt(alpha)/length(X)) # Computation of a log-convex and increasing density ratio with respect to # the Gamma distribution with shape alpha and rate beta.. # Pure computation, no graphical displays... { # If X ~ P, where P(dx) = e^theta(x) P_alpha,beta(dx) # with P_alpha,beta the distribution of a Gamma dsitribution # with shape alpha and rate beta, then beta*X ~ Pt, where # Pt(dx) = e^thetat(x) P_alpha,1(dx) and thetat(x) = theta(x/beta). # Hence we multiply the data by beta, estimate thetat and obtain # theta(x) = thetat(x*beta) in the end. X <- X*beta # Preparations of raw data: x <- unique(sort(X)) n <- length(x) W <- W/sum(W) w <- rep(0,n) for (i in 1:n) { tmp <- (X == x[i]) w[i] <- sum(W[tmp]) } # MLE for theta in V_{0} \cap \Theta: mu <- sum(w*x) kappa <- max(1 - alpha/mu,0) # \hat{\kappa}^+ ckappa <- -alpha*log(1 - kappa) # c(\hat{\kappa}^+) LL <- -ckappa + kappa*mu - exp(-ckappa)/(1-kappa)^alpha + 1 tmp <- Optimality1.2B(x,w,mu,alpha) if (length(tmp$tnew) > 0){ ht0 <- max(tmp$htnew) t0 <- tmp$tnew[which.max(tmp$htnew)] }else{ ht0 <- 0 } # Bookkeeping: nr.local.search <- 0 nr.Newton <- 0 if (ht0 < delta2){ print(paste('Parametric fit sufficient: Gamma(', alpha,',',beta,'). Log-likelihood = 0.',sep='')) return(list(tau=NULL,theta=NULL)) } tau <- 0 theta <- c(-ckappa,kappa) # Check for global optimality: tmp <- Optimality2.2B(x,w,tau,theta,alpha) if (length(tmp$tnew) > 0){ ht0 <- max(tmp$htnew) t0 <- tmp$tnew[which.max(tmp$htnew)] }else{ ht0 <- 0 } while (ht0 > delta2){ nr.local.search <- nr.local.search + 1 # Store current value of LL: LL.old <- LL # Add a new knot t0 to tau: j <- sum(tau <= t0) m <- length(tau) if (j == 0){ # Then tau[1] > 0 tau <- c(t0,tau) theta <- c(theta[1],theta) } if (j > 0 && j < m){ theta1 <- ((tau[j+1] - t0)*theta[j] + (t0 - tau[j])*theta[j+1])/(tau[j+1] - tau[j]) tau <- c(tau[1:j],t0,tau[(j+1):m]) theta <- c(theta[1:j],theta1,theta[(j+1):(m+1)]) } if (j == m){ theta1 <- theta[m] + (t0 - tau[m])*theta[m+1] tau <- c(tau,t0) theta <- c(theta[1:m],theta1,theta[m+1]) } wt <- LocalLinearSplines2.2B(tau,x,w) proposal <- LocalNewton.2B(wt,tau,theta,alpha,delta1) nr.Newton <- nr.Newton + 1 while (proposal$dirderiv > delta1){ LL.Newton <- LL theta.new <- proposal$theta.new # 2nd step size correction: corr <- LocalStepSize2.2B(tau,theta,theta.new) if (length(corr$tau) < length(tau)){ tau <- corr$tau theta <- corr$theta wt <- LocalLinearSplines2.2B(tau,x,w) }else{ theta <- theta.new } # Normalization: theta <- LocalNormalize.2B(tau,theta,alpha) LL <- sum(wt*theta) if (LL > LL.Newton){ proposal <- LocalNewton.2B(wt,tau,theta,alpha,delta1) nr.Newton <- nr.Newton + 1 }else{ proposal$dirderiv <- 0 } } # Check for global optimality: if (LL > LL.old){ tmp <- Optimality2.2B(x,w,tau,theta,alpha) if (length(tmp$tnew) > 0){ ht0 <- max(tmp$htnew) t0 <- tmp$tnew[which.max(tmp$htnew)] }else{ ht0 <- 0 } }else{ ht0 <- 0 } } # Rescale to account for beta != 1 x <- x/beta tau <- tau/beta theta[length(tau)+1] <- beta*theta[length(tau)+1] return(list(tau=tau,theta=theta,x=x,w=w, LL=LL,nr.local.search=nr.local.search,nr.Newton=nr.Newton)) } ### Auxiliary programs 4: ### Checking for knew potential knots LocalDirDeriv1.2B <- function(t,x,w,alpha=1) # This procedure computes the directional derivatives # h(t[i]) := DL(theta,V_{t[i]}) # for 1 <= i <= length(t) in case of # theta(x) = 0 . # In addition, it returns the right-sided derivatives # h'(t[i]+) . # It is only used to illustrate the methods and check the # algorithms; in the estimation procedure, only the function # Optimality1.2B() is used. { h0t <- rep(0,length(t)) h1t <- rep(0,length(t)) for (i in 1:length(t)){ JJ <- (x > t[i]) h0t[i] <- sum(w[JJ]*(x[JJ] - t[i])) - alpha*(1 - pgamma(t[i],alpha + 1,1)) + t[i]*(1 - pgamma(t[i],alpha,1)) h1t[i] <- sum(w[!JJ]) - pgamma(t[i],alpha,1) } return(cbind("h(t)"=h0t,"h'(t+)"=h1t)) } Optimality1.2B <- function(x,w,mu=sum(w*x),alpha=1) # Determine knots t0 with locally maximal directional # derivatives DL(theta,V_{t0}) in case of # theta(x) = 0 { n <- length(x) t0 <- rep(NA,n) ht0 <- rep(NA,n) t0[1] <- 0 ht0[1] <- mu - alpha Wl <- 0 for (i in 1:(n-1)){ Wl <- Wl + w[i] t0tmp <- qgamma(Wl,alpha,1) if (t0tmp > x[i] && t0tmp < x[i+1]){ t0[i+1] <- t0tmp ht0[i+1] <- sum((x[(i+1):n] - t0[i+1])*w[(i+1):n]) - alpha*(1 - pgamma(t0[i+1],alpha + 1,1)) + t0[i+1]*(1 - pgamma(t0[i+1],alpha,1)) } } return(list(tnew=t0[!is.na(t0)],htnew=ht0[!is.na(t0)])) } LocalDirDeriv2.2B <- function(t,x,w,tau,theta,alpha=1) # This procedure computes the directional derivatives # h(t[i]) := DL((tau,theta),V_{t[i],(tau,theta)}) # and the right-sided derivatives # h'(t[i]+) # for 1 <= i <= length(t). # It is only used to illustrate the methods and check the # algorithms; in the estimation procedure, only the function # Optimality2.2B() is used. { m <- length(tau) h0t <- rep(0,length(t)) h1t <- rep(0,length(t)) II <- which(t < tau[1]) for (i in II){ Gatmp <- exp(theta[1])*pgamma(t[i],alpha,1) Fhtmp <- sum(w[x <= t[i]]) JJ <- (x > t[i] & x <= tau[1]) h1t[i] <- Fhtmp - Gatmp h0t[i] <- (t[i] - tau[1])*h1t[i] + exp(theta[1])*(tau[1]*G.2B(t[i],tau[1],alpha) - alpha*G.2B(t[i],tau[1],alpha + 1)) - sum(w[JJ]*(tau[1] - x[JJ])) } if (m > 1){ for (j in 1:(m-1)){ II <- which(t >= tau[j] & t < tau[j+1]) dj <- tau[j+1] - tau[j] JJ <- (x > tau[j] & x <= tau[j+1]) Hj <- sum(w[JJ]*(tau[j+1] - x[JJ]))/dj - J10.2B(theta[j],theta[j+1],tau[j],tau[j+1],alpha) for (i in II){ thstar <- ((tau[j+1] - t[i])*theta[j] + (t[i] - tau[j])*theta[j+1])/dj JJi <- (x > tau[j] & x <= t[i]) h1t[i] <- -J00.2B(theta[j],thstar,tau[j],t[i],alpha) + sum(w[JJi]) - Hj h0t[i] <- (t[i] - tau[j])* (h1t[i] + J01.2B(theta[j],thstar,tau[j],t[i],alpha)) - sum(w[JJi]*(x[JJi] - tau[j])) } } } II <- which(t >= tau[m]) for (i in II){ thstar <- theta[m] + (t[i] - tau[m])*theta[m+1] JJ <- (x > tau[m] & x <= t[i]) K0tmp <- K0.2B(thstar,theta[m+1],t[i],alpha) Phtmp <- sum(w[x > t[i]]) h1t[i] <- K0tmp - Phtmp h0t[i] <- (t[i] - tau[m])* (h1t[i] + J01.2B(theta[m],thstar,tau[m],t[i],alpha)) - sum(w[JJ]*(x[JJ] - tau[m])) } return(cbind("h(t)"=h0t,"h'(t+)"=h1t)) } Optimality2.2B <- function(x,w,tau,theta,alpha=1) # Determine knots t0 with locally maximal directional # derivatives DL(theta,V_{t0,theta}) in case of a # non-affine current fit given by (tau,theta). # In case t[1] > 0, t0 = 0 has to be considered. { t <- unique(sort(c(x,tau))) if (t[1] > 0){ # Then tau[1] > 0, so we need to consder t0 = 0, this is done in the end. checkt0_0 <- TRUE }else{ # Otherwise tau[1] = 0, so t0 = 0 does not need to be considered. checkt0_0 <- FALSE t <- t[-1] } # first we check on the intervals (t[i],t[i+1]) nt <- length(t) tnew <- rep(NA,nt-1) htnew <- rep(NA,nt-1) m <- length(tau) II <- which(t < tau[1]) for (i in II){ Fhtmp <- sum(w[x <= t[i]]) # h'(t[i]+): Gatmp.l <- exp(theta[1])*pgamma(t[i],alpha,1) h1.l <- Fhtmp - Gatmp.l # h'(t[i+1]-): Gatmp.r <- exp(theta[1])*pgamma(t[i+1],alpha,1) h1.r <- Fhtmp - Gatmp.r if (h1.l > 0 && h1.r < 0){ ttmp <- InverseF0Theta.2B(theta[1], Fhtmp, alpha) if (ttmp > t[i] && ttmp < t[i+1]){ tnew[i] <- ttmp # thstar <- theta[2] + (ttmp - tau[1])*theta[1] JJ <- (x > ttmp & x <= tau[1]) htnew[i] <- exp(theta[1])* (tau[1]*G.2B(ttmp,tau[1],alpha) - alpha*G.2B(ttmp,tau[1],alpha + 1)) - sum(w[JJ]*(tau[1] - x[JJ])) } } } if (m > 1){ for (j in 1:(m-1)){ II <- which(t >= tau[j] & t < tau[j+1]) dj <- tau[j+1] - tau[j] JJ <- (x > tau[j] & x <= tau[j+1]) Hj <- sum(w[JJ]*(tau[j+1] - x[JJ]))/dj - J10.2B(theta[j],theta[j+1],tau[j],tau[j+1],alpha) for (i in II){ JJi <- (x > tau[j] & x <= t[i]) Wi <- sum(w[JJi]) # h'(t[i]+): thstar <- ((tau[j+1] - t[i])*theta[j] + (t[i] - tau[j])*theta[j+1])/dj h1.l <- -J00.2B(theta[j],thstar,tau[j],t[i],alpha) + Wi - Hj # h'(t[i+1]-): thstar <- ((tau[j+1] - t[i+1])*theta[j] + (t[i+1] - tau[j])*theta[j+1])/dj h1.r <- -J00.2B(theta[j],thstar,tau[j],t[i+1],alpha) + Wi - Hj if (h1.l > 0 && h1.r < 0){ stmp <- (theta[j+1] - theta[j])/dj ttmp <- InverseJ00.2B(theta[j],stmp,tau[j],Wi-Hj,alpha) if (ttmp > t[i] && ttmp < t[i+1]){ tnew[i] <- ttmp thstar <- ((tau[j+1] - ttmp)*theta[j] + (ttmp - tau[j])*theta[j+1])/dj htnew[i] <- (ttmp - tau[j])* J01.2B(theta[j],thstar,tau[j],ttmp,alpha) - sum(w[JJi]*(x[JJi] - tau[j])) } } } } } II <- which(t[1:(nt-1)] >= tau[m]) for (i in II){ Shtmp <- sum(w[x > t[i]]) # h'(t[i]+): thstar.l <- theta[m] + (t[i] - tau[m])*theta[m+1] K0tmp.l <- K0.2B(thstar.l,theta[m+1],t[i],alpha) h1.l <- K0tmp.l - Shtmp # h'(t[i+1]-): thstar.r <- theta[m] + (t[i+1] - tau[m])*theta[m+1] K0tmp.r <- K0.2B(thstar.r,theta[m+1],t[i+1],alpha) h1.r <- K0tmp.r - Shtmp if (h1.l > 0 && h1.r < 0){ ttmp <- InverseK0.2B(theta[m],theta[m+1],tau[m],Shtmp,alpha) if (ttmp > t[i] && ttmp < t[i+1]){ tnew[i] <- ttmp thstar <- theta[m] + (ttmp - tau[m])*theta[m+1] JJ <- (x > tau[m] & x <= ttmp) htnew[i] <- (ttmp - tau[m])* J01.2B(theta[m],thstar,tau[m],ttmp,alpha) - sum(w[JJ]*(x[JJ] - tau[m])) } } } if (checkt0_0 == TRUE) { tnew <- c(0,tnew) JJ <- (x > 0 & x <= tau[1]) httmp <- exp(theta[1])* (tau[1]*pgamma(tau[1],alpha,1) - alpha*pgamma(tau[1],alpha + 1,1)) - sum(w[JJ]*(tau[1] - x[JJ])) htnew <- c(httmp,htnew) } return(list(tnew=tnew[!is.na(tnew)],htnew=htnew[!is.na(tnew)])) } ### Auxiliary programs 3: ### 2nd step size correction LocalConvexity.2B <- function(tau,theta) # Changes of slope, # each slope and change of slope should be > 0 and # if theta defines a strictly convex and increasing function on {tau[i] : ...}. { m <- length(tau) if (m >= 2){ dtau <- tau[2:m] - tau[1:(m-1)] dtheta <- theta[2:m] - theta[1:(m-1)] dtdt <- c(0,dtheta/dtau,theta[m+1]) } else{ dtdt <- c(0,theta[2]) } chofslope <- dtdt[2:(m+1)] - dtdt[1:m] return(chofslope) } LocalStepSize2.2B <- function(tau,theta,theta.new) # Replace theta with # (1 - t0)*theta + t0*theta.new , # where t0 is the largest number in [0,1] such that # the new theta defines still a convex function. # Return new tuples tau and theta, possibly shorter # than the original ones. { m <- length(tau) chofsl <- LocalConvexity.2B(tau,theta) chofsl.new <- LocalConvexity.2B(tau,theta.new) if (m == 1 && chofsl.new <= 0){ theta <- c(0,0) # Constant function equal to 0 with an active knot at tau[1] return(list(tau=tau,theta=theta)) } if (min(chofsl.new) <= 0){ tmp <- which(chofsl.new <= 0) t0 <- chofsl[tmp]/(chofsl[tmp] - chofsl.new[tmp]) j0 <- which.min(t0) t0 <- t0[j0] j0 <- tmp[j0] theta <- (1 - t0)*theta + t0*theta.new tau <- tau[-j0] theta <- theta[-j0] m <- m-1 chofsl <- LocalConvexity.2B(tau,theta) while (m > 1 && min(chofsl) <= 0){ j0 <- which.min(chofsl) tau <- tau[-j0] theta <- theta[-j0] m <- m-1 chofsl <- LocalConvexity.2B(tau,theta) } if (m == 1 && chofsl <= 0){ theta <- c(0,0) # Constant function equal to 0 with an active knot at tau[1] } }else{ theta <- theta.new } return(list(tau=tau,theta=theta)) } ### Auxiliary programs 2: ### Normalization, log-likelihood plus derivatives ### and Newton step with 1st step size correction LocalNormalize.2B <- function(tau,theta,alpha=1) # Normalize theta such that it defines # a log-probability density { m <- length(tau) integral <- exp(theta[1])*pgamma(tau[1], alpha, 1) + K0.2B(theta[m],theta[m+1],tau[m],alpha) if (m > 1) { integral <- integral + sum(J00.2B(theta[1:(m-1)],theta[2:m], tau[1:(m-1)],tau[2:m],alpha)) } theta.new <- theta theta.new[1:m] <- theta.new[1:m] - log(integral) return(theta.new) } LocalLL.2B <- function(wt,tau,theta,alpha=1) # Log-likelihood { m <- length(tau) integral <- exp(theta[1])*pgamma(tau[1], alpha, 1) + K0.2B(theta[m],theta[m+1],tau[m],alpha) if (m > 1) { integral <- integral + sum(J00.2B(theta[1:(m-1)],theta[2:m], tau[1:(m-1)],tau[2:m],alpha)) } LL <- sum(wt*theta) + 1 - integral return(LL) } LocalLL2.2B <- function(wt,tau,theta,alpha=1) # Log-likelihood plus its gradient vector # and minus Hessian matrix { m <- length(tau) integral <- exp(theta[1])*pgamma(tau[1], alpha, 1) + K0.2B(theta[m],theta[m+1],tau[m],alpha) if (m > 1) { integral <- integral + sum(J00.2B(theta[1:(m-1)],theta[2:m],tau[1:(m-1)],tau[2:m],alpha)) } LL <- sum(wt*theta) + 1 - integral GLL <- wt GLL[1] <- GLL[1] - exp(theta[1])*pgamma(tau[1], alpha, 1) GLL[m] <- GLL[m] - K0.2B(theta[m],theta[m+1],tau[m],alpha) GLL[m+1] <- GLL[m+1] - K1.2B(theta[m],theta[m+1],tau[m],alpha) if (m > 1) { GLL[1:(m-1)] <- GLL[1:(m-1)] - J10.2B(theta[1:(m-1)],theta[2:m],tau[1:(m-1)],tau[2:m],alpha) GLL[2:m] <- GLL[2:m] - J01.2B(theta[1:(m-1)],theta[2:m],tau[1:(m-1)],tau[2:m],alpha) } MHLL <- matrix(0,m+1,m+1) MHLL[1,1] <- exp(theta[1])*pgamma(tau[1], alpha, 1) MHLL[m,m] <- MHLL[m,m] + K0.2B(theta[m],theta[m+1],tau[m],alpha) MHLL[m,m+1] <- K1.2B(theta[m],theta[m+1],tau[m],alpha) MHLL[m+1,m] <- MHLL[m,m+1] MHLL[m+1,m+1] <- K2.2B(theta[m],theta[m+1],tau[m],alpha) if (m > 1){ A20 <- cbind(1:(m-1),1:(m-1)) MHLL[A20] <- MHLL[A20] + J20.2B(theta[1:(m-1)],theta[2:m], tau[1:(m-1)],tau[2:m],alpha) A02 <- cbind(2:m,2:m) MHLL[A02] <- MHLL[A02] + J02.2B(theta[1:(m-1)],theta[2:m], tau[1:(m-1)],tau[2:m],alpha) tmp11 <- J11.2B(theta[1:(m-1)],theta[2:m], tau[1:(m-1)],tau[2:m],alpha) MHLL[cbind(1:(m-1),2:m)] <- tmp11 MHLL[cbind(2:m,1:(m-1))] <- tmp11 } return(list(LL=LL,GLL=GLL,MHLL=MHLL)) } LocalNewton.2B <- function(wt,tau,theta,alpha=1,delta0=10^(-11)) # Newton step with (1st) step size correction { m <- length(tau) integral <- exp(theta[1])*pgamma(tau[1], alpha, 1) + K0.2B(theta[m],theta[m+1],tau[m],alpha) if (m > 1) { integral <- integral + sum(J00.2B(theta[1:(m-1)],theta[2:m],tau[1:(m-1)],tau[2:m],alpha)) } LL <- sum(wt*theta) + 1 - integral GLL <- wt GLL[1] <- GLL[1] - exp(theta[1])*pgamma(tau[1], alpha, 1) GLL[m] <- GLL[m] - K0.2B(theta[m],theta[m+1],tau[m],alpha) GLL[m+1] <- GLL[m+1] - K1.2B(theta[m],theta[m+1],tau[m],alpha) if (m > 1) { GLL[1:(m-1)] <- GLL[1:(m-1)] - J10.2B(theta[1:(m-1)],theta[2:m],tau[1:(m-1)],tau[2:m],alpha) GLL[2:m] <- GLL[2:m] - J01.2B(theta[1:(m-1)],theta[2:m],tau[1:(m-1)],tau[2:m],alpha) } MHLL <- matrix(0,m+1,m+1) MHLL[1,1] <- exp(theta[1])*pgamma(tau[1], alpha, 1) MHLL[m,m] <- MHLL[m,m] + K0.2B(theta[m],theta[m+1],tau[m],alpha) MHLL[m,m+1] <- K1.2B(theta[m],theta[m+1],tau[m],alpha) MHLL[m+1,m] <- MHLL[m,m+1] MHLL[m+1,m+1] <- K2.2B(theta[m],theta[m+1],tau[m],alpha) if (m > 1){ A20 <- cbind(1:(m-1),1:(m-1)) MHLL[A20] <- MHLL[A20] + J20.2B(theta[1:(m-1)],theta[2:m], tau[1:(m-1)],tau[2:m],alpha) A02 <- cbind(2:m,2:m) MHLL[A02] <- MHLL[A02] + J02.2B(theta[1:(m-1)],theta[2:m], tau[1:(m-1)],tau[2:m],alpha) tmp11 <- J11.2B(theta[1:(m-1)],theta[2:m], tau[1:(m-1)],tau[2:m],alpha) MHLL[cbind(1:(m-1),2:m)] <- tmp11 MHLL[cbind(2:m,1:(m-1))] <- tmp11 } dtheta <- qr.solve(MHLL,GLL) dirderiv <- sum(dtheta * GLL) theta.new <- theta + dtheta LL.new <- LocalLL.2B(wt,tau,theta.new,alpha) while (LL.new < LL + dirderiv/3 && dirderiv > delta0){ theta.new <- (theta + theta.new)/2 dirderiv <- dirderiv/2 LL.new <- LocalLL.2B(wt,tau,theta.new,alpha) } return(list(theta.new=theta.new,dirderiv=dirderiv, LL=LL,LL.new=LL.new)) } ### Auxiliary programs 1: ### Linear splines and ### Basic functions J(.,.), K(.,.) plus 1st and 2nd ### partial derivatives LocalLinearSplines1.2B <- function(tau,x=tau) # For a vector tau with m >= 1 strictly increasing # components, this function computes a design matrix Dx # for the set of piecewise linear functions f on the real # line with kinks only on {tau[j] : 1 <= j <= m}, evaluated # at x. The basis corresponds to the parameter vector theta # with components # theta[j] = f(tau[j]), 1 <= j <= m, # theta[m+1] = f'(tau[m] +) . { m <- length(tau) n <- length(x) Dx <- matrix(0,n,m+1) tmp <- (x <= tau[1]) Dx[tmp,1] <- 1 if (m > 1){ dtau <- tau[2:m] - tau[1:(m-1)] for (j in 1:(m-1)){ tmp <- (x > tau[j] & x <= tau[j+1]) Dx[tmp,j ] <- (tau[j+1] - x[tmp])/dtau[j] Dx[tmp,j+1] <- (x[tmp] - tau[j])/dtau[j] } } tmp <- (x > tau[m]) Dx[tmp,m] <- 1 Dx[tmp,m+1] <- x[tmp] - tau[m] return(Dx) } LocalLinearSplines2.2B <- function(tau, x=tau,w=rep(1/length(x),length(x))) # For linear splines with given knots # tau[1] < tau[2] < ... < tau[m] # and arbitrary vectors x and w of equal length with w >= 0, # this procedure returns a weight vector wt, such that for # the basis functions f_1, f_2, ..., f_m, f_{m+1}, f_{m+2} # described/computed in LocalLinearSplines1(), # wt[j] = sum_{i=1}^n w[i]*f_j(x[i]) . { Dx <- LocalLinearSplines1.2B(tau,x) wt <- colSums(Dx * w) return(wt) } G.2B <- function(a,b,s=1) # G(a,b;s) = F0(b,s) - F0(a,s). { pgamma(b, s, 1) - pgamma(a, s, 1) } J00.2B <- function(x,y,a=0,b=1,alpha=1,epsilon=1e-4) # J00.2B(x,y,a,b) # = integral_a^b # exp((1-v(z))x + v(z)y)*pgamma(z, alpha, 1) dz # with v(z) = (z - a)/(b - a). { xt <- (b*x - a*y)/(b - a + (a == b)) yt <- (y - x)/(b - a + (a == b)) bt <- (1 - yt)*b at <- (1 - yt)*a J00xy <- rep(0,length(x)) II <- (yt < 1) J00xy[II] <- exp(xt[II] + log(G.2B(at[II], bt[II], alpha))) * (1 - yt[II])^(- alpha) II <- (yt >= 1) # In this case, use an upper bound (exact value in case of yt = 1) J00xy[II] <- exp(xt[II] + (yt[II] - 1)*b[II]) * (b[II]^alpha - a[II]^alpha) / gamma(alpha + 1) return(J00xy) # # An alternative option for numerical stability. # # If it is used, the last "return(J00xy)" should be deleted. # II <- (yt <= 1-epsilon | yt >= 1) # if (all(II)) { # return(J00xy) # }else{ # II <- !II # JJ <- (1:length(xt))[II] # c0 <- (b[II]^ alpha - a[II]^ alpha) / alpha # c1 <- (b[II]^(alpha + 1) - a[II]^(alpha + 1))/(alpha + 1) # J00.L <- exp(xt[II] + (yt[II] - 1)*c1/c0) * # c0/gamma(alpha) # J00.U <- exp(xt[II] + (yt[II] - 1)*a[II]) * # c0/gamma(alpha) # KK <- which(J00xy[II] < J00.L | J00.U < J00xy[II]) # JJ <- JJ[KK] # if (length(JJ) > 0){ # J00xy[JJ] <- (J00.L[KK] + J00.U[KK])/2 # } # return(J00xy) # } } J10.2B <- function(x,y,a=0,b=1,alpha=1) # J10.2B(x,y,a,b) # = integral_a^b # (1-v(z))*exp((1-v(z))x + v(z)y)*pgamma(z, alpha, 1) dz # with v(z) = (z - a)/(b - a). { xt <- (b*x - a*y)/(b - a + (a == b)) yt <- (y - x)/(b - a + (a == b)) bt <- (1 - yt)*b at <- (1 - yt)*a J10xy <- exp(xt) * (1 - yt)^(- alpha - 1) * (bt * G.2B(at, bt, alpha) - alpha*G.2B(at, bt, alpha + 1)) / (b - a + (a == b)) return(J10xy) } J01.2B <- function(x,y,a=0,b=1,alpha=1) # J01.2B(x,y,a,b) # = integral_a^b # v(z)*exp((1-v(z))x + v(z)y)*pgamma(z, alpha, 1) dz # with v(z) = (z - a)/(b - a). { xt <- (b*x - a*y)/(b - a + (a == b)) yt <- (y - x)/(b - a + (a == b)) bt <- (1 - yt)*b at <- (1 - yt)*a J01xy <- exp(xt) * (1 - yt)^(- alpha - 1) * (-at * G.2B(at, bt, alpha) + alpha*G.2B(at, bt, alpha + 1)) / (b - a + (a == b)) return(J01xy) } J20.2B <- function(x,y,a=0,b=1,alpha=1) # J20.2B(x,y,a,b) # = integral_a^b # (1-v(z))^2*exp((1-v(z))x + v(z)y)*pgamma(z, alpha, 1) dz # with v(z) = (z - a)/(b - a). { xt <- (b*x - a*y)/(b - a + (a == b)) yt <- (y - x)/(b - a + (a == b)) bt <- (1 - yt)*b at <- (1 - yt)*a J20xy <- exp(xt) * (1 - yt)^(- alpha - 2) * ( bt^2 * G.2B(at, bt, alpha ) -2*alpha*bt * G.2B(at, bt, alpha + 1) +alpha*(alpha+1) * G.2B(at, bt, alpha + 2)) / (b - a + (a == b))^2 return(J20xy) } J11.2B <- function(x,y,a=0,b=1,alpha=1) # J11.2B(x,y,a,b) # = integral_a^b # (1-v(z))*v(z)*exp((1-v(z))x + v(z)y)*pgamma(z, alpha, 1) dz # with v(z) = (z - a)/(b - a). { xt <- (b*x - a*y)/(b - a + (a == b)) yt <- (y - x)/(b - a + (a == b)) bt <- (1 - yt)*b at <- (1 - yt)*a J11xy <- exp(xt) * (1 - yt)^(- alpha - 2) * ( -at*bt * G.2B(at, bt, alpha ) +alpha *(at+bt) * G.2B(at, bt, alpha + 1) -alpha*(alpha+1) * G.2B(at, bt, alpha + 2)) / (b - a + (a == b))^2 return(J11xy) } J02.2B <- function(x,y,a=0,b=1,alpha=1) # J02.2B(x,y,a,b) # = integral_a^b # v(z)^2*exp((1-v(z))x + v(z)y)*pgamma(z, alpha, 1) dz # with v(z) = (z - a)/(b - a). { xt <- (b*x - a*y)/(b - a + (a == b)) yt <- (y - x)/(b - a + (a == b)) bt <- (1 - yt)*b at <- (1 - yt)*a J02xy <- exp(xt) * (1 - yt)^(- alpha - 2) * ( at^2 * G.2B(at, bt, alpha ) - 2*alpha*at * G.2B(at, bt, alpha + 1) + alpha*(alpha+1) * G.2B(at, bt, alpha + 2)) / (b - a + (a == b))^2 return(J02xy) } K0.2B <- function(x,y,a=0,alpha=1) # K0.2B(x,y,a) = integral_a^Inf exp(x + y*(z-a))*pgamma(z, alpha, 1) dz { if(y >= 1) { return(Inf) } return(exp(x - a*y)*(1-pgamma((1-y)*a, alpha, 1))/(1-y)^alpha) } K1.2B <- function(x,y,a=0,alpha=1) # K1.2B(x,y,a) = dK0.2B(x,y,a)/dy { at <- (1 - y)*a return(exp(x - a*y)* ( -at*(1 - pgamma(at, alpha, 1)) +alpha*(1 - pgamma(at, alpha + 1, 1)))/ (1 - y)^(alpha + 1)) } K2.2B <- function(x,y,a=0,alpha=1) # K2.2B(x,y,a) = d^2K0.2B(x,y,a)/dy^2 { at <- (1 - y)*a return(exp(x - a*y)* ( at^2*(1 - pgamma(at, alpha , 1)) -2*alpha*at*(1 - pgamma(at, alpha + 1, 1)) +alpha*(alpha+1)*(1 - pgamma(at, alpha + 2, 1)))/ (1 - y)^(alpha + 2)) } InverseF0Theta.2B <- function(theta0, c, alpha=1) # Determines the unique solution tau of # exp(theta0)*F0(tau,alpha) = c. { d <- c*exp(-theta0) if (d < 0 || d >= 1) { return(NULL) } return(qgamma(d, alpha, 1)) } InverseJ00.2B <- function(theta0, theta1, tau0, c, alpha=1) # Determines the unique solution tau of # J00(theta0, theta0 + theta1*(tau-tau0),tau0,tau) = c. { if (theta1 < 0 || theta1 >= 1) { return(NULL) } d <- c*(1-theta1)^alpha*exp(theta1*tau0 - theta0) + pgamma((1 - theta1)*tau0, alpha, 1) if (d >= 1 || d < 0) { return(NULL) } return(qgamma(d, alpha, 1)/(1 - theta1)) } InverseK0.2B <- function(theta0, theta1, tau0, c, alpha=1) # Determines the unique solution tau of # K0.2B(theta0 + theta1*(tau-tau0),theta1,tau) = c. { if (theta1 < 0 || theta1 >= 1) { return(NULL) } d <- c*(1-theta1)^alpha*exp(theta1*tau0 - theta0) if (d > 1 || d <= 0) { return(NULL) } return(qgamma(1-d, alpha, 1)/(1 - theta1)) } LocalInterpolate <- function(x,d) # For a vector x with n = length(x) >= 2, this function returns # a vector xx with n*d - d + 1 components, such that # xx[(j*d-d+1):(j*d)] == x[j] + (x[j+1] - x[j])*(1:d)/d # for 1 <= j < n. { n <- length(x) xx <- rep(NA,n*d-d+1) xx[1] <- x[1] for (k in 1:d){ lambda <- k/d xx[1 + k + (0:(n-2))*d] <- (1 - lambda)*x[1:(n-1)] + lambda*x[2:n] } return(xx) } ### Auxiliary programs 0: ### Data simulation and ### plotting functions Evaluate.2B <- function(t, tau, theta) # For plotting purposes only. # For a vector t of nt = length(t) non-negative elements, # this function returns the value of theta(t). An active # knot at 0 is artificially added in case tau[1] > 0. { if (tau[1] > 0) { tau <- c(0,tau) theta <- c(theta[1], theta) } m <- length(tau) dtau <- tau[2:m] - tau[1:(m-1)] dtheta <- theta[2:m] - theta[1:(m-1)] dtdt <- c(dtheta/dtau,theta[m+1]) tau <- c(tau, Inf) nt <- length(t) res <- rep(0, nt) for(i in 1:nt) { j <- which(tau[1:m] <= t[i] & t[i] < tau[2:(m+1)]) res[i] <- theta[j] + (t[i]-tau[j])*dtdt[j] } return(res) } Simulate.2B <- function(n,tau,theta,alpha=1,beta=1,sim.num=FALSE) # The function returns simulated data # (X_1, X_2, ..., X_N) in R^N with density # f_theta(x) = f_o(x) * exp(theta(x)) / M(theta) # where f_o(x) is the density of Gamma(alpha,beta) and # theta is a piecewise linear convex and increasing # function, with values theta at knots tau. The last # element of theta is the last slope of theta(x) and it has # to be smaller than beta. # If sim.num is TRUE, returns the number of needed # simulations and the vectors tau and theta. { # Add a knot at 0 if needed if(tau[1] > 0){ tau <- c(0,tau) theta <- c(theta[1], theta) } m <- length(tau) if(is.infinite(tau[m])) { tau <- tau[1:(m-1)] m <- m - 1 } # Knots, slopes, values of theta at the knots: if (m == 1) { theta_prime <- theta[m+1] } else { theta_prime <- rep(0, m) theta_prime[1:(m-1)] <- (theta[2:m]-theta[1:(m-1)])/(tau[2:m]-tau[1:(m-1)]) theta_prime[m] <- theta[m+1] } tau <- c(tau, Inf) theta <- theta[1:m] beta_b <- beta - theta_prime[m] # g function, g(x) in [0,1] for all x g <- function(x){ k_x <- (1:m)[(tau[1:m] <= x & x < tau[2:(m+1)])] res <- (theta_prime[k_x] - theta_prime[m])*(x - tau[k_x]) - theta_prime[m]*tau[k_x] + theta[k_x] - theta[1] return(exp(res)) } # Acceptance rejection: X <- rep(0,n) i <- 1 nn <- 0 while(i <= n){ Y <- rgamma(1, shape = alpha, rate = beta_b) U <- runif(1) if(U <= g(Y)){ X[i] <- Y i <- i+1 } nn <- nn + 1 } if(sim.num == TRUE) return(list(X = X,sim.num = nn,theta = c(theta, theta_prime[m]),tau = tau)) else return(X) }