### Main programs:

TailInflation.A0 <- function(X,W=rep(1/length(X),length(X)),
	xmin=1.1*min(X)-0.1*max(X),
	xmax=1.1*max(X)-0.1*min(X),
	delta1=10^(-10)/length(X),
	delta2=10^(-3)/length(X))
# Computation of a log-convex density ratio with respect to
# the standard Gaussian distribution.
# Graphical display of all intermediate steps...
{
	# 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])
	}
	# Gaussian MLE theta (for variance 1):
	mu <- sum(w*x)
	LL <- sum(w*(mu*x - mu^2/2))
	tmp <- Optimality1.2A(x,w,mu)
	if (length(tmp$tnew) > 0){
		ht0 <- max(tmp$htnew)
		t0 <- tmp$tnew[which.max(tmp$htnew)]
	}else{
		ht0 <- 0
	}
		
	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(mu*c(xmin,xmax) - mu^2/2),
		main=paste('Starting point: ',
			'LL = ',round(LL,digits=9),sep=''))
	lines(c(xmin,xmax),mu*c(xmin,xmax)-mu^2/2,lwd=2)
	rug(x)

	# Bookkeeping:
	nr.local.search <- 0
	nr.Newton <- 0
	
	if (ht0 < delta2){
		print(paste('Parametric fit sufficient: N(',
			round(mu,4),', 1)',sep=''))
		return(list(tau=NULL,theta=NULL,mu=mu))
	}

	# 1st local search:
	nr.local.search <- nr.local.search + 1
	tau <- t0
	theta <- c(mu,mu*t0 - mu^2/2,mu)
	wt <- LocalLinearSplines2.2A(tau,x,w)
	proposal <- LocalNewton.2A(wt,tau,theta,delta1)
	nr.Newton <- nr.Newton + 1
	# For graphical display:
	xtx <- c(xmin,tau,xmax)
	Dtmp <- LocalLinearSplines1.2A(tau,xtx)
	while (proposal$dirderiv > delta1){
		theta.new <- proposal$theta.new
		# Graphical display:
		message <- readline("New proposal:")
		thetaxtx     <- Dtmp %*% theta
		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(c(thetaxtx,thetaxtx.new)),
			main=paste('LL =',
				round(LL,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')
		rug(x)
		# Update of theta:
		theta <- theta.new
		theta <- LocalNormalize.2A(tau,theta)
		proposal <- LocalNewton.2A(wt,tau,theta,delta1)
		nr.Newton <- nr.Newton + 1
		LL <- proposal$LL
	}
	# 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)))
	points(tau,thetaxtx[2:(length(tau)+1)])
	lines(xtx,thetaxtx,lwd=2)
	rug(x)
	
	# Check for global optimality:
	tmp <- Optimality2.2A(x,w,tau,theta)
	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){
			theta1 <- theta[2] + (t0 - tau[1])*theta[1]
			tau <- c(t0,tau)
			theta <- c(theta[1],theta1,theta[2:(m+2)])
		}
		if (j > 0 && j < m){
			theta1 <- ((tau[j+1] - t0)*theta[j+1] +
				(t0 - tau[j])*theta[j+2])/(tau[j+1] - tau[j])
			tau <- c(tau[1:j],t0,tau[(j+1):m])
			theta <- c(theta[1:(j+1)],theta1,theta[(j+2):(m+2)])
		}
		if (j == m){
			theta1 <- theta[m+1] + (t0 - tau[m])*theta[m+2]
			tau <- c(tau,t0)
			theta <- c(theta[1:(m+1)],theta1,theta[m+2])
		}
		# Graphical display:
		message <- readline("New knot:")
		xtx <- c(xmin,tau,xmax)
		Dtmp <- LocalLinearSplines1.2A(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=t0,col='blue')
		rug(x)

		wt <- LocalLinearSplines2.2A(tau,x,w)
		proposal <- LocalNewton.2A(wt,tau,theta,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=t0,col='blue')
			rug(x)
			# 2nd step size correction:
			corr <- LocalStepSize2.2A(tau,theta,theta.new)
			if (length(corr$tau) < length(tau)){
				tau <- corr$tau
				theta <- corr$theta
				wt <- LocalLinearSplines2.2A(tau,x,w)
				# Graphical display:
				message <- readline("2nd step size correction:")
				xtx <- c(xmin,tau,xmax)
				Dtmp <- LocalLinearSplines1.2A(tau,xtx)
				thetaxtx <- Dtmp %*% theta
				lines(xtx,thetaxtx,lwd=2,col='magenta')
			}else{
				theta <- theta.new
			}
			# Normalization:
			theta <- LocalNormalize.2A(tau,theta)
			LL <- sum(wt*theta)
			if (LL > LL.Newton){
				proposal <- LocalNewton.2A(wt,tau,theta,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)])
		rug(x)
		
		# Check for global optimality:
		if (LL > LL.old){
			tmp <- Optimality2.2A(x,w,tau,theta)
			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)])
	rug(x)
	
	return(list(tau=tau,theta=theta,x=x,w=w,
		LL=LL,nr.local.search=nr.local.search,nr.Newton=nr.Newton))
}


TailInflation.A <- function(X,W=rep(1/length(X),length(X)),
	delta1=10^(-10)/length(X),
	delta2=10^(-3)/length(X))
# Computation of a log-convex density ratio with respect to
# the standard Gaussian distribution.
# Pure computation, no graphical displays...
{
	# 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])
	}
	# Gaussian MLE theta (for variance 1):
	mu <- sum(w*x)
	LL <- sum(w*(mu*x - mu^2/2))
	tmp <- Optimality1.2A(x,w,mu)
	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: N(',
			round(mu,4),', 1)',sep=''))
		return(list(tau=NULL,theta=NULL,mu=mu))
	}

	# 1st local search:
	nr.local.search <- nr.local.search + 1
	tau <- t0
	theta <- c(mu,mu*t0 - mu^2/2,mu)
	wt <- LocalLinearSplines2.2A(tau,x,w)
	proposal <- LocalNewton.2A(wt,tau,theta,delta1)
	nr.Newton <- nr.Newton + 1
	while (proposal$dirderiv > delta1){
		theta.new <- proposal$theta.new
		# Update of theta:
		theta <- theta.new
		theta <- LocalNormalize.2A(tau,theta)
		proposal <- LocalNewton.2A(wt,tau,theta,delta1)
		nr.Newton <- nr.Newton + 1
		LL <- proposal$LL
	}
	
	# Check for global optimality:
	tmp <- Optimality2.2A(x,w,tau,theta)
	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){
			theta1 <- theta[2] + (t0 - tau[1])*theta[1]
			tau <- c(t0,tau)
			theta <- c(theta[1],theta1,theta[2:(m+2)])
		}
		if (j > 0 && j < m){
			theta1 <- ((tau[j+1] - t0)*theta[j+1] +
				(t0 - tau[j])*theta[j+2])/(tau[j+1] - tau[j])
			tau <- c(tau[1:j],t0,tau[(j+1):m])
			theta <- c(theta[1:(j+1)],theta1,theta[(j+2):(m+2)])
		}
		if (j == m){
			theta1 <- theta[m+1] + (t0 - tau[m])*theta[m+2]
			tau <- c(tau,t0)
			theta <- c(theta[1:(m+1)],theta1,theta[m+2])
		}

		wt <- LocalLinearSplines2.2A(tau,x,w)
		proposal <- LocalNewton.2A(wt,tau,theta,delta1)
		nr.Newton <- nr.Newton + 1
		while (proposal$dirderiv > delta1){
			LL.Newton <- LL
			theta.new <- proposal$theta.new
			# 2nd step size correction:
			corr <- LocalStepSize2.2A(tau,theta,theta.new)
			if (length(corr$tau) < length(tau)){
				tau <- corr$tau
				theta <- corr$theta
				wt <- LocalLinearSplines2.2A(tau,x,w)
			}else{
				theta <- theta.new
			}
			# Normalization:
			theta <- LocalNormalize.2A(tau,theta)
			LL <- sum(wt*theta)
			if (LL > LL.Newton){
				proposal <- LocalNewton.2A(wt,tau,theta,delta1)
				nr.Newton <- nr.Newton + 1
			}else{
				proposal$dirderiv <- 0	
			}
		}
		
		# Check for global optimality:
		if (LL > LL.old){
			tmp <- Optimality2.2A(x,w,tau,theta)
			if (length(tmp$tnew) > 0){
				ht0 <- max(tmp$htnew)
				t0 <- tmp$tnew[which.max(tmp$htnew)]
			}else{
				ht0 <- 0
			}
		}else{
			ht0 <- 0
		}
	}
		
	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.2A <- function(t,x,w,mu=sum(x*w))
# This procedure computes the directional derivatives
#    h(t[i]) := DL(theta,V_{t[i]})
# for 1 <= i <= length(t) in case of
#    theta(x) = mu*x - mu^2/2 .
# 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.2A() 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])) -
			K1.2A(mu*t[i] - mu^2/2,mu,t[i])
		h1t[i] <- sum(w[!JJ]) - pnorm(t[i] - mu)
	}
	return(cbind("h(t)"=h0t,"h'(t+)"=h1t))
}

Optimality1.2A <- function(x,w,mu=sum(w*x))
# Determine knots t0 with locally maximal directional
# derivatives DL(theta,V_{t0}) in case of
# theta(x) = mu*x - mu^2/2.
{
	n <- length(x)
	t0 <- rep(NA,n-1)
	ht0 <- rep(NA,n-1)
	Wl <- 0
	for (i in 1:(n-1)){
		Wl <- Wl + w[i]
		t0tmp <- mu + qnorm(Wl)
		if (t0tmp > x[i] && t0tmp < x[i+1]){
			t0[i] <- t0tmp
			ht0[i] <- sum((x[(i+1):n] - t0[i])*w[(i+1):n]) -
				K1.2A(mu*t0[i] - mu^2/2,mu,t0[i])
		}
	}
	return(list(tnew=t0[!is.na(t0)],htnew=ht0[!is.na(t0)]))
}

LocalDirDeriv2.2A <- function(t,x,w,tau,theta)
# 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.2A() is used.
{
	m <- length(tau)
	h0t <- rep(0,length(t))
	h1t <- rep(0,length(t))
	II <- which(t < tau[1])
	for (i in II){
		thstar <- theta[2] + (t[i] - tau[1])*theta[1]
		K0tmp <- K0.2A(thstar,-theta[1],-t[i])
		Fhtmp <- sum(w[x <= t[i]])
		JJ <- (x > t[i] & x <= tau[1])
		h1t[i] <- Fhtmp - K0tmp
		h0t[i] <- (t[i] - tau[1])*
			(h1t[i] - J10.2A(thstar,theta[2],t[i],tau[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.2A(theta[j+1],theta[j+2],tau[j],tau[j+1])
			for (i in II){
				thstar <- ((tau[j+1] - t[i])*theta[j+1] +
					(t[i] - tau[j])*theta[j+2])/dj
				JJi <- (x > tau[j] & x <= t[i])
				h1t[i] <- -J00.2A(theta[j+1],thstar,tau[j],t[i]) +
					sum(w[JJi]) - Hj
				h0t[i] <- (t[i] - tau[j])*
					(h1t[i] + J01.2A(theta[j+1],thstar,tau[j],t[i])) -
					sum(w[JJi]*(x[JJi] - tau[j]))
			}
		}
	}
	II <- which(t >= tau[m])
	for (i in II){
		thstar <- theta[m+1] + (t[i] - tau[m])*theta[m+2]
		JJ <- (x > tau[m] & x <= t[i])
		K0tmp <- K0.2A(thstar,theta[m+2],t[i])
		Phtmp <- sum(w[x > t[i]])
		h1t[i] <- K0tmp - Phtmp
		h0t[i] <- (t[i] - tau[m])*
			(h1t[i] + J01.2A(theta[m+1],thstar,tau[m],t[i])) -
			sum(w[JJ]*(x[JJ] - tau[m]))
	}
	return(cbind("h(t)"=h0t,"h'(t+)"=h1t))
}


Optimality2.2A <- function(x,w,tau,theta)
# 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).
{
	t <- unique(sort(c(x,tau)))
	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]+):
		thstar.l <- theta[2] + (t[i] - tau[1])*theta[1]
		K0tmp.l <- K0.2A(thstar.l,-theta[1],-t[i])
		h1.l <- Fhtmp - K0tmp.l
		# h'(t[i+1]-):
		thstar.r <- theta[2] + (t[i+1] - tau[1])*theta[1]
		K0tmp.r <- K0.2A(thstar.r,-theta[1],-t[i+1])
		h1.r <- Fhtmp - K0tmp.r
		if (h1.l > 0 && h1.r < 0){
			ttmp <- InverseK0.2A(theta[2],theta[1],tau[1],Fhtmp,-1)
			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] <- (tau[1] - ttmp)*
						J10.2A(thstar,theta[2],ttmp,tau[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.2A(theta[j+1],theta[j+2],tau[j],tau[j+1])
			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+1] +
					(t[i] - tau[j])*theta[j+2])/dj
				h1.l <- -J00.2A(theta[j+1],thstar,tau[j],t[i]) +
					Wi - Hj
				# h'(t[i+1]-):
				thstar <- ((tau[j+1] - t[i+1])*theta[j+1] +
					(t[i+1] - tau[j])*theta[j+2])/dj
				h1.r <- -J00.2A(theta[j+1],thstar,tau[j],t[i+1]) +
					Wi - Hj
				if (h1.l > 0 && h1.r < 0){
					stmp <- (theta[j+2] - theta[j+1])/dj
					ttmp <- InverseJ00.2A(theta[j+1],stmp,tau[j],Wi-Hj)
					if (ttmp > t[i] && ttmp < t[i+1]){
						tnew[i] <- ttmp
						thstar <- ((tau[j+1] - ttmp)*theta[j+1] +
							(ttmp - tau[j])*theta[j+2])/dj
						htnew[i] <- (ttmp - tau[j])*	
							J01.2A(theta[j+1],thstar,tau[j],ttmp) -
							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+1] + (t[i] - tau[m])*theta[m+2]
		K0tmp.l <- K0.2A(thstar.l,theta[m+2],t[i])
		h1.l <- K0tmp.l - Shtmp
		# h'(t[i+1]-):
		thstar.r <- theta[m+1] + (t[i+1] - tau[m])*theta[m+2]
		K0tmp.r <- K0.2A(thstar.r,theta[m+2],t[i+1])
		h1.r <- K0tmp.r - Shtmp
		if (h1.l > 0 && h1.r < 0){
			ttmp <- InverseK0.2A(theta[m+1],theta[m+2],tau[m],Shtmp,+1)
			if (ttmp > t[i] && ttmp < t[i+1]){
				tnew[i] <- ttmp
				thstar <- theta[m+1] + (ttmp - tau[m])*theta[m+2]
				JJ <- (x > tau[m] & x <= ttmp)
				htnew[i] <- (ttmp - tau[m])*
						J01.2A(theta[m+1],thstar,tau[m],ttmp) -
					sum(w[JJ]*(x[JJ] - tau[m]))
			}
		}		
	}
	return(list(tnew=tnew[!is.na(tnew)],htnew=htnew[!is.na(tnew)]))
}

### Auxiliary programs 3:
### 2nd step size correction

LocalConvexity.2A <- function(tau,theta)
# Changes of slope,
# all of which should be > 0 if theta defines
# a strictly concave function on {tau[i] : ...}.
{
	m <- length(tau)
	if (m >= 2){
		dtau <- tau[2:m] - tau[1:(m-1)]
		dtheta <- theta[3:(m+1)] - theta[2:m]
		dtdt <- c(theta[1],dtheta/dtau,theta[m+2])
	}
	else{
		dtdt <- c(theta[1],theta[3])
	}
	chofslope <- dtdt[2:(m+1)] - dtdt[1:m]
	return(chofslope)
}

LocalStepSize2.2A <- 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.2A(tau,theta)
	chofsl.new <- LocalConvexity.2A(tau,theta.new)
	if (m == 1 && chofsl.new <= 0){
		t0 <- chofsl/(chofsl - chofsl.new)
		theta <- (1 - t0)*theta + t0*theta.new
		theta[c(1,3)] <- mean(theta[c(1,3)])
		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+1)]
		m <- m-1
		chofsl <- LocalConvexity.2A(tau,theta)
		while (m > 1 && min(chofsl) <= 0){
			j0 <- which.min(chofsl)
			tau <- tau[-j0]
			theta <- theta[-(j0+1)]
			m <- m-1
			chofsl <- LocalConvexity.2A(tau,theta)
		}
		if (m == 1 && chofsl <= 0){
			theta[c(1,3)] <- mean(theta[c(1,3)])
		}
	}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.2A <- function(tau,theta)
# Normalize theta such that it defines
# a log-probability density
{
	m <- length(tau)
	integral <- K0.2A(theta[2],-theta[1],-tau[1]) +
		K0.2A(theta[m+1], theta[m+2], tau[m])
	if (m > 1)
	{
		integral <- integral +
			sum(J00.2A(theta[2:m],theta[3:(m+1)],
				tau[1:(m-1)],tau[2:m]))
	}
	theta.new <- theta
	theta.new[2:(m+1)] <- theta.new[2:(m+1)] - log(integral)
	return(theta.new)
}

LocalLL.2A <- function(wt,tau,theta)
# Log-likelihood
{
	m <- length(tau)
	integral <- K0.2A(theta[2],-theta[1],-tau[1]) +
		K0.2A(theta[m+1], theta[m+2], tau[m])
	if (m > 1)
	{
		integral <- integral +
			sum(J00.2A(theta[2:m],theta[3:(m+1)],
				tau[1:(m-1)],tau[2:m]))
	}
	LL <- sum(wt*theta) + 1 - integral
	return(LL)
}

LocalLL2.2A <- function(wt,tau,theta)
# Log-likelihood plus its gradient vector
# and minus Hessian matrix
{
	m <- length(tau)
	integral <- K0.2A(theta[2],-theta[1],-tau[1]) +
		K0.2A(theta[m+1], theta[m+2], tau[m])
	if (m > 1){
		integral <- integral +
			sum(J00.2A(theta[2:m],theta[3:(m+1)],
				tau[1:(m-1)],tau[2:m]))
	}
	LL <- sum(wt*theta) + 1 - integral
	GLL <- wt
	GLL[1] <- GLL[1] + K1.2A(theta[2], -theta[1], -tau[1])
	GLL[2] <- GLL[2] - K0.2A(theta[2], -theta[1], -tau[1])
	GLL[m+1] <- GLL[m+1] - K0.2A(theta[m+1], theta[m+2], tau[m])
	GLL[m+2] <- GLL[m+2] - K1.2A(theta[m+1], theta[m+2], tau[m])
	if (m > 1){
		GLL[2:m] <- GLL[2:m] -
			J10.2A(theta[2:m],theta[3:(m+1)],tau[1:(m-1)],tau[2:m])
		GLL[3:(m+1)] <- GLL[3:(m+1)] -
			J01.2A(theta[2:m],theta[3:(m+1)],tau[1:(m-1)],tau[2:m])
	}
	MHLL <- matrix(0,m+2,m+2)
	MHLL[1,1] <- K2.2A(theta[2], - theta[1], -tau[1])
	MHLL[1,2] <- - K1.2A(theta[2], - theta[1], -tau[1])
	MHLL[2,1] <- MHLL[1,2]
	MHLL[2,2] <- K0.2A(theta[2], -theta[1], -tau[1])
	MHLL[m+1,m+1] <- MHLL[m+1,m+1] +
		K0.2A(theta[m+1], theta[m+2], tau[m])
	MHLL[m+1,m+2] <- K1.2A(theta[m+1], theta[m+2], tau[m])
	MHLL[m+2,m+1] <- MHLL[m+1,m+2]
	MHLL[m+2,m+2] <- K2.2A(theta[m+1], theta[m+2], tau[m])
	if (m > 1){
		A20 <- cbind(2:m,2:m)
		MHLL[A20] <- MHLL[A20] +
			J20.2A(theta[2:m],theta[3:(m+1)],
				tau[1:(m-1)],tau[2:m])
		A02 <- cbind(3:(m+1),3:(m+1))
		MHLL[A02] <- MHLL[A02] +
			J02.2A(theta[2:m],theta[3:(m+1)],
				tau[1:(m-1)],tau[2:m])
		tmp11 <- J11.2A(theta[2:m],theta[3:(m+1)],
			tau[1:(m-1)],tau[2:m])
		MHLL[cbind(2:m,3:(m+1))] <- tmp11
		MHLL[cbind(3:(m+1),2:m)] <- tmp11
	}
	return(list(LL=LL,GLL=GLL,MHLL=MHLL))
}

LocalNewton.2A <- function(wt,tau,theta,
	delta0=10^(-11))
# Newton step with (1st) step size correction
{
	m <- length(tau)
	integral <- K0.2A(theta[2],-theta[1],-tau[1]) +
		K0.2A(theta[m+1], theta[m+2], tau[m])
	if (m > 1){
		integral <- integral +
			sum(J00.2A(theta[2:m],theta[3:(m+1)],
				tau[1:(m-1)],tau[2:m]))
	}
	LL <- sum(wt*theta) + 1 - integral
	GLL <- wt
	GLL[1] <- GLL[1] + K1.2A(theta[2], -theta[1], -tau[1])
	GLL[2] <- GLL[2] - K0.2A(theta[2], -theta[1], -tau[1])
	GLL[m+1] <- GLL[m+1] - K0.2A(theta[m+1], theta[m+2], tau[m])
	GLL[m+2] <- GLL[m+2] - K1.2A(theta[m+1], theta[m+2], tau[m])
	if (m > 1){
		GLL[2:m] <- GLL[2:m] -
			J10.2A(theta[2:m],theta[3:(m+1)],tau[1:(m-1)],tau[2:m])
		GLL[3:(m+1)] <- GLL[3:(m+1)] -
			J01.2A(theta[2:m],theta[3:(m+1)],tau[1:(m-1)],tau[2:m])
	}
	MHLL <- matrix(0,m+2,m+2)
	MHLL[1,1] <- K2.2A(theta[2], - theta[1], -tau[1])
	MHLL[1,2] <- - K1.2A(theta[2], - theta[1], -tau[1])
	MHLL[2,1] <- MHLL[1,2]
	MHLL[2,2] <- K0.2A(theta[2], -theta[1], -tau[1])
	MHLL[m+1,m+1] <- MHLL[m+1,m+1] +
		K0.2A(theta[m+1], theta[m+2], tau[m])
	MHLL[m+1,m+2] <- K1.2A(theta[m+1], theta[m+2], tau[m])
	MHLL[m+2,m+1] <- MHLL[m+1,m+2]
	MHLL[m+2,m+2] <- K2.2A(theta[m+1], theta[m+2], tau[m])
	if (m > 1){
		A20 <- cbind(2:m,2:m)
		MHLL[A20] <- MHLL[A20] +
			J20.2A(theta[2:m],theta[3:(m+1)],
				tau[1:(m-1)],tau[2:m])
		A02 <- cbind(3:(m+1),3:(m+1))
		MHLL[A02] <- MHLL[A02] +
			J02.2A(theta[2:m],theta[3:(m+1)],
				tau[1:(m-1)],tau[2:m])
		tmp11 <- J11.2A(theta[2:m],theta[3:(m+1)],
			tau[1:(m-1)],tau[2:m])
		MHLL[cbind(2:m,3:(m+1))] <- tmp11
		MHLL[cbind(3:(m+1),2:m)] <- tmp11
	}
	dtheta <- qr.solve(MHLL,GLL)
	dirderiv <- sum(dtheta * GLL)
	theta.new <- theta + dtheta
	LL.new <- LocalLL.2A(wt,tau,theta.new)
	while (LL.new < LL + dirderiv/3 && dirderiv > delta0){
		theta.new <- (theta + theta.new)/2
		dirderiv <- dirderiv/2
		LL.new <- LocalLL.2A(wt,tau,theta.new)
	}
	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.2A <- 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[1]   = f'(tau[1] -) ,
#    theta[j+1] = f(tau[j]), 1 <= j <= m ,
#    theta[m+2] = f'(tau[m] +) .
{
	m <- length(tau)
	n <- length(x)
	Dx <- matrix(0,n,m+2)
	tmp <- (x <= tau[1])
	Dx[tmp,1] <- x[tmp] - tau[1]
	Dx[tmp,2] <- 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+1] <- (tau[j+1] - x[tmp])/dtau[j]
			Dx[tmp,j+2] <- (x[tmp] - tau[j])/dtau[j]
		}
	}
	tmp <- (x > tau[m])
	Dx[tmp,m+1] <- 1
	Dx[tmp,m+2] <- x[tmp] - tau[m]
	return(Dx)
}


LocalLinearSplines2.2A <- 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.2A(),
#    wt[j] = sum_{i=1}^n w[i]*f_j(x[i]) .
{
	Dx <- LocalLinearSplines1.2A(tau,x)
	wt <- colSums(Dx * w)
	return(wt)
}


J00.2A <- function(x,y,a=0,b=1)
# J00.2A(x,y,a,b)
#    =  integral_a^b
#          exp((1-v(z))x + v(z)y)*dnorm(z) dz
# with v(z) = (z - a)/(b - a).
{
	xt <- (b*x - a*y)/(b - a + (a == b))
	yt <- (y - x)/(b - a + (a == b))
	at <- a - yt
	bt <- b - yt
	mt <- (a + b)/2 - yt
	dt <- (b - a)/2
	pdiff <- pnorm(bt) - pnorm(at)
	tmp <- (at > 0)
	pdiff[tmp] <- pnorm(-at[tmp]) - pnorm(-bt[tmp])
	J00xy <- exp(xt + yt^2/2 +
		pmax(log(pdiff),-mt^2/2 + log(2*pnorm(dt) - 1)))
	return(J00xy)
}

J10.2A <- function(x,y,a=0,b=1)
# J10.2A(x,y,a,b)
#    =  integral_a^b
#          (1-v(z))*exp((1-v(z))x + v(z)y)*dnorm(z) dz
# with v(z) = (z - a)/(b - a).
{
	xt <- (b*x - a*y)/(b - a + (a == b))
	yt <- (y - x)/(b - a + (a == b))
	at <- a - yt
	bt <- b - yt
	pdiff <- pnorm(bt) - pnorm(at)
	tmp <- (at > 0)
	pdiff[tmp] <- pnorm(-at[tmp]) - pnorm(-bt[tmp])
	J10xy <- exp(xt + yt^2/2)*
		(bt*pdiff + dnorm(bt) - dnorm(at))/
			(b - a + (a == b))
	return(J10xy)
}

J01.2A <- function(x,y,a=0,b=1)
# J01.2A(x,y,a,b)
#    =  integral_a^b
#          v(z)*exp((1-v(z))x + v(z)y)*dnorm(z) dz
# with v(z) = (z - a)/(b - a).
# Symmetry considerations reveal that
#    J01.2A(x,y,a,b) = J10.2A(y,x,-b,-a)
{
	return(J10.2A(y,x,-b,-a))
}

J20.2A <- function(x,y,a=0,b=1)
# J20.2A(x,y,a,b)
#    =  integral_a^b
#          (1-v(z))^2*exp((1-v(z))x + v(z)y)*dnorm(z) dz
# with v(z) = (z - a)/(b - a).
{
	xt <- (b*x - a*y)/(b - a + (a == b))
	yt <- (y - x)/(b - a + (a == b))
	at <- a - yt
	bt <- b - yt
	pdiff <- pnorm(bt) - pnorm(at)
	tmp <- (at > 0)
	pdiff[tmp] <- pnorm(-at[tmp]) - pnorm(-bt[tmp])
	J20xy <- exp(xt + yt^2/2)*
		((1 + bt^2)*pdiff +
		 (at - 2*bt)*dnorm(at) + bt*dnorm(bt))/
			(b - a + (a == b))^2
	return(J20xy)
}

J11.2A <- function(x,y,a=0,b=1)
# J11.2A(x,y,a,b)
#    =  integral_a^b
#          (1-v(z))*v(z)*exp((1-v(z))x + v(z)y)*dnorm(z) dz
# with v(z) = (z - a)/(b - a).
{
	xt <- (b*x - a*y)/(b - a + (a == b))
	yt <- (y - x)/(b - a + (a == b))
	at <- a - yt
	bt <- b - yt
	pdiff <- pnorm(bt) - pnorm(at)
	tmp <- (at > 0)
	pdiff[tmp] <- pnorm(-at[tmp]) - pnorm(-bt[tmp])
	J11xy <- exp(xt + yt^2/2)*
		(-(1 + at*bt)*pdiff + bt*dnorm(at) - at*dnorm(bt))/
			(b - a + (a == b))^2
	return(J11xy)
}

J02.2A <- function(x,y,a=0,b=1)
# J02.2A(x,y,a,b)
#    =  integral_a^b
#          v(z)^2*exp((1-v(z))x + v(z)y)*dnorm(z) dz
# with v(z) = (z - a)/(b - a).
{
	xt <- (b*x - a*y)/(b - a + (a == b))
	yt <- (y - x)/(b - a + (a == b))
	at <- a - yt
	bt <- b - yt
	pdiff <- pnorm(bt) - pnorm(at)
	tmp <- (at > 0)
	pdiff[tmp] <- pnorm(-at[tmp]) - pnorm(-bt[tmp])
	J02xy <- exp(xt + yt^2/2)*
		((1 + at^2)*pdiff +
		 (2*at - bt)*dnorm(bt) - at*dnorm(at))/
			(b - a + (a == b))^2
	return(J02xy)
}


K0.2A <- function(x,y,a=0)
# K0.2A(x,y,a) = integral_a^Inf exp(x + y*(z-a))*dnorm(z) dz
{
	return(exp(x - a*y + y^2/2)*pnorm(y-a))
}

K1.2A <- function(x,y,a=0)
# K1.2A(x,y,a) = dK0.2A(x,y,a)/dy
{
	return(exp(x - a*y + y^2/2)*
		((y - a)*pnorm(y-a) + dnorm(y - a)))
}

K2.2A <- function(x,y,a=0)
# K2.2A(x,y,a) = d^2K0.2A(x,y,a)/dy^2
{
	return(exp(x - a*y + y^2/2)*
		((1 + (y - a)^2)*pnorm(y-a) +
			(y - a)*dnorm(y - a)))
}


InverseK0.2A <- function(theta0, theta1, tau0, c, s=1)
# Determines the unique solution tau of
#    K0.2A(theta0 + theta1*(tau-tau0),s*theta1,s*tau) = c.
{
	if (c <= 0)
	{
		return(NULL)
	}
	d <- c*exp(-theta0 + theta1*tau0 - theta1^2/2)
	if (d >= 1)
	{
		return(NULL)
	}
	return(theta1 - s*qnorm(d))
}

InverseJ00.2A <- function(theta0, theta1, tau0, c)
# Determines the unique solution tau of
#    J00.2A(theta0, theta0 + theta1*(tau-tau0),tau0,tau) = c.
{
	d <- pnorm(tau0 - theta1) +
		c*exp(-theta0 + theta1*tau0 - theta1^2/2)
	if (d >= 1 || d <= 0)
	{
		return(NULL)
	}
	return(theta1 + qnorm(d))
}


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)
}