R: Calculate the log likelihood of a smoothing spline given the...
likelihood.smooth.spline
R Documentation
Calculate the log likelihood of a smoothing spline given the data
Description
Calculate the (log) likelihood of a spline given the data used to fit
the spline, g. The likelihood consists of two main parts:
1) (weighted) residuals sum of squares, and 2) a penalty term. The
penalty term consists of a smoothing parameterlambda
and a roughness measure of the spline
J(g) = int g''(t) dt. Hence, the overall log likelihood is
log L(g|x) = (y-g(x))'W(y-g(x)) + λ J(g)
In addition to the overall likelihood, all its seperate
components are also returned.
Note: when fitting a smooth spline with (x,y) values where the
x's are not unique, smooth.spline will replace
such (x,y)'s with a new pair (x,y') where y' is a
reweighted average on the original y's. It is important to
be aware of this. In such cases, the resulting smooth.spline
object does not contain all (x,y)'s and therefore this
function will not calculate the weighted residuals sum of square on
the original data set, but on the data set with unique x's.
See examples below how to calculate the likelihood for the spline with
the original data.
Usage
## S3 method for class 'smooth.spline'
likelihood(object, x=NULL, y=NULL, w=NULL, base=exp(1),
rel.tol=.Machine$double.eps^(1/8), ...)
Arguments
object
The smooth.spline object.
x, y
The x and y values for which the (weighted) likelihood will
be calculated. If x is of type xy.coords any value of
argument y will be omitted. If x==NULL, the x and y values
of the smoothing spline will be used.
w
The weights for which the (weighted) likelihood will be
calculated. If NULL, weights equal to one are assumed.
base
The base of the logarithm of the likelihood. If NULL,
the non-logged likelihood is returned.
rel.tol
The relative tolerance used in the call to
integrate.
...
Not used.
Details
The roughness penalty for the smoothing spline, g, fitted
from data in the interval [a,b] is defined as
J(g) = int_a^b g''(t) dt
which is the same as
J(g) = g'(b) - g'(a)
The latter is calculated internally by using
predict.smooth.spline.
Value
Returns the overall (log) likelihood of class
SmoothSplineLikelihood, a class with the following attributes:
wrss
the (weighted) residual sum of square
penalty
the penalty which is equal to -lambda*roughness.
lambda
the smoothing parameter
roughness
the value of the roughness functional given the
specific smoothing spline and the range of data
Author(s)
Henrik Bengtsson
See Also
smooth.spline and robustSmoothSpline().
Examples
# Define f(x)
f <- expression(0.1*x^4 + 1*x^3 + 2*x^2 + x + 10*sin(2*x))
# Simulate data from this function in the range [a,b]
a <- -2; b <- 5
x <- seq(a, b, length=3000)
y <- eval(f)
# Add some noise to the data
y <- y + rnorm(length(y), 0, 10)
# Plot the function and its second derivative
plot(x,y, type="l", lwd=4)
# Fit a cubic smoothing spline and plot it
g <- smooth.spline(x,y, df=16)
lines(g, col="yellow", lwd=2, lty=2)
# Calculating the (log) likelihood of the fitted spline
l <- likelihood(g)
cat("Log likelihood with unique x values:\n")
print(l)
# Note that this is not the same as the log likelihood of the
# data on the fitted spline iff the x values are non-unique
x[1:5] <- x[1] # Non-unique x values
g <- smooth.spline(x,y, df=16)
l <- likelihood(g)
cat("\nLog likelihood of the *spline* data set:\n");
print(l)
# In cases with non unique x values one has to proceed as
# below if one want to get the log likelihood for the original
# data.
l <- likelihood(g, x=x, y=y)
cat("\nLog likelihood of the *original* data set:\n");
print(l)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(aroma.light)
aroma.light v3.2.0 (2016-01-06) successfully loaded. See ?aroma.light for help.
> png(filename="/home/ddbj/snapshot/RGM3/R_BC/result/aroma.light/likelihood.smooth.spline.Rd_%03d_medium.png", width=480, height=480)
> ### Name: likelihood.smooth.spline
> ### Title: Calculate the log likelihood of a smoothing spline given the
> ### data
> ### Aliases: likelihood.smooth.spline
> ### Keywords: methods smooth internal
>
> ### ** Examples
>
> # Define f(x)
> f <- expression(0.1*x^4 + 1*x^3 + 2*x^2 + x + 10*sin(2*x))
>
> # Simulate data from this function in the range [a,b]
> a <- -2; b <- 5
> x <- seq(a, b, length=3000)
> y <- eval(f)
>
> # Add some noise to the data
> y <- y + rnorm(length(y), 0, 10)
>
> # Plot the function and its second derivative
> plot(x,y, type="l", lwd=4)
>
> # Fit a cubic smoothing spline and plot it
> g <- smooth.spline(x,y, df=16)
> lines(g, col="yellow", lwd=2, lty=2)
>
> # Calculating the (log) likelihood of the fitted spline
> l <- likelihood(g)
>
> cat("Log likelihood with unique x values:\n")
Log likelihood with unique x values:
> print(l)
Likelihood of smoothing spline: -293865.1
Log base: 2.718282
Weighted residuals sum of square: 293865.2
Penalty: -0.1210083
Smoothing parameter lambda: 0.0009257147
Roughness score: 130.7187
>
> # Note that this is not the same as the log likelihood of the
> # data on the fitted spline iff the x values are non-unique
> x[1:5] <- x[1] # Non-unique x values
> g <- smooth.spline(x,y, df=16)
> l <- likelihood(g)
>
> cat("\nLog likelihood of the *spline* data set:\n");
Log likelihood of the *spline* data set:
> print(l)
Likelihood of smoothing spline: -293228.5
Log base: 2.718282
Weighted residuals sum of square: 293228.6
Penalty: -0.1210751
Smoothing parameter lambda: 0.0009261969
Roughness score: 130.7229
>
> # In cases with non unique x values one has to proceed as
> # below if one want to get the log likelihood for the original
> # data.
> l <- likelihood(g, x=x, y=y)
> cat("\nLog likelihood of the *original* data set:\n");
Log likelihood of the *original* data set:
> print(l)
Likelihood of smoothing spline: -293864.2
Log base: 2.718282
Weighted residuals sum of square: 293864.3
Penalty: -0.1210758
Smoothing parameter lambda: 0.0009261969
Roughness score: 130.7236
>
>
>
>
>
>
>
>
>
>
> dev.off()
null device
1
>