
x1 <- seq(0, 1, 0.00001)
	
fold.val <- function(x, mu, vr, thresh, smooth) {

	a <- pnorm(x1, mean=thresh, sd=smooth)
	d <- dnorm(x1, mean=mu, sd=sqrt(vr))
	y <- d*(1-a)

	yr <- rev(d*a)
	xr <- rev(2*thresh-x1)

	yr1 <- approx(xr, yr, x1)$y

	yr1[is.na(yr1)] <- 0
	return(approx(x1, y+yr1, x)$y)
}

fold.err <- function(p, xx, yy) {

	r <- fold.val(xx, mu=p[1], vr=p[2], thresh=p[3], smooth=p[4])
	return(sum((r-yy)^2))
}


h <- read.table("lnsrch.dat", col.names=c('x', 'y'))

param <- c(.413, 0.00687, .5228, .01255)
up    <- c(.45,  0.0072,  .53,   .014)
lo    <- c(.4,   0.0060,  .52,   .011)

m <- optim(param, fold.err, xx=h$x, yy=h$y, method='L-BFGS-B', lower=lo, upper=up)

cat(m$message, "\n")

plot(h$x, h$y, type='l', xlab='Distance', ylab='Density')
lines(h$x, fold.val(h$x, mu=m$par[1], vr=m$par[2], thresh=m$par[3], smooth=m$par[4]), lty=3, lwd=2, col='orange')
