oldpar <- par(no.readonly=TRUE)
# evaluation points
x = seq(-2,2,len=11)
# evaluate a standard normal distribution function
p = pnorm(x)
# combine with 1-p
mnormp = cbind(p,1-p)
# convert to surprisal values
mnorms = -log2(mnormp)
# plot the surprisal values
matplot(x, mnorms, type="l", lty=c(1,1), col=c(1,1),
ylab="Surprisal (2-bits)")
# add some log-normal error
mnormdata = exp(log(mnorms) + rnorm(11)*0.1)
# set up a b-spline basis object
nbasis = 7
sbasis = create.bspline.basis(c(-2,2),nbasis)
# define an initial coefficient matrix
cmat = matrix(0,7,1)
# set up a fdPar object for suprisal smoothing
sfd = fd(cmat, sbasis)
sfdPar = fdPar(sfd, Lfd=2, lambda=0)
# smooth the noisy data
result = smooth.surp(x, mnormdata, cmat, sfdPar)
# plot the data and the fits of the two surprisal curves
xfine = seq(-2,2,len=51)
sfine = eval.surp(xfine, result$Wfd)
matplot(xfine, sfine, type="l", lty=c(1,1), col=c(1,1))
points(x, mnormdata[,1])
points(x, mnormdata[,2])
# convert the surprisal fit values to probabilities
pfine = 2^(-sfine)
# check that they sum to one
apply(pfine,1,sum)
par(oldpar)
Run the code above in your browser using DataLab