just to clarify how I see the error, it was the mis-definition of the
penalty term in the function dv. The following code corrects this error.
What is actually being minimised at this step is the penalised deviance
conditional on the smoothing parameter. A second issue is that the optim
default ("Nelder-Mead") would have to be recycled several times to
approximate the minimum obtained by gam, however, in this example, "BFGS"
gets it straight away.
sorry for all the confusion!
greg
set.seed(1)
library(mgcv)
x<-runif(100)
lp<-exp(-2*x)*sin(8*x)
y<-rpois(100,exp(lp))
m1<-gam(y~s(x),poisson)
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)
s<-m1$sp
dv<-function(t)
{
f<-exp(X%*%t)
-2*sum(y*log(f)-f-ifelse(y==0,0,y*log(y))+y)+s*t%*%S%*%t
}
dv(b)
m1$dev+m1$sp*b%*%S%*%b
t1<-optim(rep(0,10),dv,method="BFGS")
t1$p
b
dv(t1$p)
dv(b)
fit1<-exp(X%*%t1$p)
plot(x,y)
points(x,exp(lp),pch=16,col="green3")
points(x,fitted(m1),pch=16,cex=0.5,col="blue")
points(x,fit1,cex=1.5,col="red")
> 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.
>>
>>
>
>
______________________________________________
R-help@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.
