Re: [R] gam and optim
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.
Re: [R] gam and optim
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.
[R] gam and optim
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.
