[R] Curve fitting, probably splines

Simon Wood s.wood at bath.ac.uk
Mon Apr 16 09:33:04 CEST 2012


You could use the summation convention built into mgcv:gam for this. See 
?linear.functional.terms
for details, but here is some example code, both for the exact match, 
you describe, and a noisy version. best, Simon

library(mgcv)
f2 <- function(x) 0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10 ## test func
x <- (1:1000-.5)/1000  ## quadrature points
f <- f2(x) ## function values at x
X <- matrix(x,20,1000,byrow=TRUE)
L <- X*0
for (i in 1:20) L[i,1:50+(i-1)*50] <- 1/50
y <- L%*%f ## average function values in 20 intervals

## estimate f from average values, using summation convention
## in gam. Basically smooth is evaluated at each value in X
## to give a matrix s(X) and then rowSums(s(X)*L) is predictor
## for y ...

b <- gam(y~s(X,by=L,fx=TRUE,k=20))
plot(b,shift=coef(b)[1]) ## estimate
lines(x,f,col=2)         ## truth

## noisy version...
y <- y + rnorm(10)*.1 ## add noise to averages
b <- gam(y~s(X,by=L),method="REML")
plot(b,shift=coef(b)[1]) ## estimate
lines(x,f,col=2)              ## truth



On 04/12/2012 02:45 PM, Michael Haenlein wrote:
> Dear all,
>
> This is probably more related to statistics than to [R] but I hope someone
> can give me an idea how to solve it nevertheless:
>
> Assume I have a variable y that is a function of x: y=f(x). I know the
> average value of y for different intervals of x. For example, I know that
> in the interval[0;x1] the average y is y1, in the interval [x1;x2] the
> average y is y2 and so forth.
>
> I would like to find a line of minimum curvature so that the average values
> of y in each interval correspond to y1, y2, ...
>
> My idea was to use (cubic) splines. But the problem I have seems somewhat
> different to what is usually done with splines. As far as I understand it,
> splines help to find a curve that passes a set of given points. But I don't
> have any points, I only have average values of y per interval.
>
> If you have any suggestions on how to solve this, I'd love to hear them.
>
> Thanks very much in advance,
>
> Michael
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



More information about the R-help mailing list