Positive smoothing with the fda-package (Functional data analysis)

Question

In the book Functional Data Analysis with R (Ramsay&Silverman) there is described the possibility to do the "positive smoothing" if it’s needed instead of the "normal smoothing".

In the books Chapter 5.4.1 ist says:

...function w(t) is now the logarithm of the data-fitting function x(t) = exp[w(t)], and consequently is unconstrained as to its sign, while at the same time the fitting function is guaranteed to be positive. It can go as close to zero as we like by permitting the values of w(t) to be arbitrarily large negative numbers. For example, we can smooth Vancouver’s mean daily precipitation data, which can have zero but not negative values, using these commands using the function smooth.pos in R.

So as I understand, if I do smoothing with the function smooth.pos the smoothed function would not be negative or have negative values? I computed the example in the book with the Vancouvers Precipitation and the positive smoothed function is sometimes negative. (Picture smooth.pos)

"normal smooth" here I plotted the fd object.

"positive smooth" here I plotted the Wfd object.

Here the code to try it yourself.

### Constrained Smoothing from fda
##positive smoothing
#
library(fda)

VancPrec = CanadianWeather$dailyAv[,'Vancouver','Precipitation.mm']
time = 1:365

sum(VancPrec==0)

# saturated Fourier basis
fbasis = create.fourier.basis(c(0,365),365)

# harmonic acceleration penalty
harmaccelLfd = vec2Lfd(c(0,2*pi/365,0),c(0,365))
vanprecpar = fdPar(fbasis,harmaccelLfd,lambda=1e3)

# direct smoothing 
normalsmooth = smooth.basis(time,VancPrec,vanprecpar)
plot(normalsmooth$fd, ylim=c(-1,7)) #(graphic: Unrestricted smooth)
abline(h=0,lty = 2, col=4)

# positive smoothing
??smooth.pos

positivesmooth = smooth.pos(time,VancPrec,vanprecpar)
plot(positivesmooth$Wfdobj) #(graphic: Positive smooth)

##################plot after evaluating

#compare the smooths
plot(time,VancPrec,col=2,xlab='day',ylab='precipitation',
     cex.lab=1.5,cex.axis=1.5)
lines(normalsmooth$fd,col=4,lwd=2)

# There is no direct plot command for positive fd objects, but the
# eval.posfd function will do the exponentiation for you. 
??eval.fd
presvals = eval.fd(time,normalsmooth$fd)
plot(presvals, ylim=c(-1,7))
abline(h=0,lty = 2, col=4)

??eval.posfd
posvals = eval.posfd(time,positivesmooth$Wfdobj)
plot(posvals, ylim=c(-1,7))
abline(h=0,lty = 2, col=4)
#lines(time,posvals,col=5,lwd=2)

So why is it then negetaive? Also I tried the smooth.pos on my data and it produced some functions completely in the negative region. Thx.

Answer 1

The simple answer to the "why is it then negative" is that $\lambda$ is not large enough. If you increase the $\lambda$ from 1e3 to 1e4 you will get a correct answer.

That being said, even the documentation example does not provide positive smoothed estimates if we use lambda = 1e3 instead of 1e4. Because of this I am uncertain that smooth.pos actually is working correctly; this behaviour looks like a bug to me. I would not use that function at this point. Maybe consider using some other packages from the Functional Data Analysis CRAN Task View.

If we need strictly non-negative smoothed data I would suggest treating the data as coming from a distribution with a non-negative domain to begin with. Then we can, for example, use mgcv::gam, treat our data as coming from a Gamma distribution and get a non-negative smooth mean function in a straightforward and coherent way.

Answer 2

Look at the function not only being unexpectedly negative, but also the scale is entirely different. The reason is that the output is the logarithm.

So you need to convert it back by exponentiating.

y_out = exp(eval.fd(time,positivesmooth$Wfdobj))
plot(time, y_out, type="l", ylim = c(-1,max(y_out)), ylab="value")
lines(time, rep(0, length(time)), lty=2,col=4)
points(time, VancPrec, pch=21, col=1, bg=1, cex=0.3)   

Answer 3

If you are interested in positive smoothing, another solution could be to re-parameterize the coefficients of your smoother (I think it is what the fda pkg does as well). A similar idea is clearly illustrated in Unimodal smoothing (see log-concave smoothing - section 4 - in appendix you will also find some Matlab codes).

So to solve our problem we need two ingredients:

• The smoother - I think we can use a Whittaker smoother (see e.g. A Perfect Smoother). In this case my basis is just the identity matrix
• The parameterization - I will use $\exp$ parameterization

What follows can be easily adapted if you want to use P-splines (see Matlab codes in the paper I mentioned before) and other penalized smoothing methods.

Our smoothing problem can then be written as $$ \min_{z} S = \|y - \mu \|^{2} + \lambda \|D^{d}z\|^2 $$ where $\mu = \exp(z)$, $z$ is a vector of unknown and $D^{d}$ is a difference operator of order $d$ which is used as penalty to regularize the behavior of adjacent $z$s (let say $d$ = 3). Note that this parameterization ensures positive estimates whatever value of $\lambda$. The $\min$ problem can be solved in many ways (Eilers' paper uses IWLS).

The results for your example with a $\lambda = 5e4$ should look like in the plot below (which gives a reasonably smooth fit within the settings above)