[R] gam and optim

[Greg Dropkin]

please ignore this, I see the error.

greg

> hi
>
> probably a silly mistake, but I expected gam to minimise the penalised
> deviance.
>
> thanks
>
> greg
>
> set.seed(1)
> library(mgcv)
> x<-runif(100)
> lp<-exp(-2*x)*sin(8*x)
> y<-rpois(100,exp(lp))
> plot(x,y)
> m1<-gam(y~s(x),poisson)
> points(x,exp(lp),pch=16,col="green3")
> points(x,fitted(m1),pch=16,cex=0.5,col="blue")
> W<-diag(fitted(m1))
> X<-predict(m1,type="lpmatrix")
> S<-m1$smooth[[1]]$S[[1]]
> S<-rbind(0,cbind(0,S))
> A<-X%*%solve(t(X)%*%W%*%X+m1$sp*S)%*%t(X)%*%W
> sum(diag(A))
> sum(m1$edf)
> fit<-fitted(m1)
> b<-m1$coef
> range(exp(X%*%b)-fit)
> z<-y/fit-1+X%*%b
> range(A%*%z-X%*%b)
>
> dv<-function(t)
> {
> f<-exp(X%*%t)
> -2*sum(y*log(f)-f-ifelse(y==0,0,y*log(y))+y)+t%*%S%*%t
> }
> dv(b)
> m1$dev+b%*%S%*%b
>
> #so far, so good
>
>
> t1<-optim(rep(0,10),dv)
> t1$p
> b
>
> #different
>
> dv(t1$p)
> dv(b)
>
> #different, and dv(t1$p) is lower!
>
> fit1<-exp(X%*%t1$p)
> points(x,fit1,pch=16,cex=0.5,col="red")
>
> # different
> # gam found b which does approximate the true curve, but does not minimise
> the penalised deviance, by a long shot.
>
>