[R] How to get the penalized log likelihood from smooth.spline()?

Henrik Bengtsson hb at maths.lth.se
Wed Feb 20 17:49:46 CET 2002


I use smooth.spline(x, y) in package modreg and I would like to get
value of penalized log likelihood and preferable also its two parts. To
make clear what I am asking for (and make sure that I am asking for the
right thing) I clarify my problem trying to use the same notation as in
help(smooth.spline):

I want to find the natural cubic spline f(x) such that

   L(f) = \sum_{k=1}{n} w[k](y[k] - f(x[k])^2 + \lambda J(f)

is minimized, where J(f) := \int f''(t)^2 dt is the quadratic roughness
functional use. Since J(f) is quadratic one can find a matrix \Sigma such
that J(g) = c^T{\Sigma}c where c is the vector of spline coefficients.
With J(f) defined as above the elements of \Sigma becomes

  \Sigma_{ij} = \int \beta_i''(t)\beta_j''(t) dt

where \beta(t) is the vector of B-spline base functions. Finally, writing
the matrix W as W := diag(\sqrt{w}) one can write L(f) as

  L(f) = (y - f)^T W^2 (y - f) + \lambda c^T{\Sigma}c

which is the form used in help(smooth.spline). So back to my question,
using smooth.spline(), how can I get 

  1) the value of L(f(x)) for a my fitted (x,y) data,
  2) the value of roughness penalty c^T{\Sigma}c,
  3) the value of (y - f)^T W^2 (y - f)?

Does smooth.spline() return any of these directly or indirectly?

Thanks

Henrik Bengtsson

Dept. of Mathematical Statistics @ Centre for Mathematical Sciences 
Lund Institute of Technology/Lund University, Sweden (+2h UTC)
Office: P316, +46 46 222 9611 (phone), +46 46 222 4623 (fax)
h b @ m a t h s . l t h . s e
http://www.maths.lth.se/bioinformatics/





-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._



More information about the R-help mailing list