Cross validation from Jacobian matrix of likelihood for glm(poisson, identity)

The difference of Standard Error between glm(y~x, family=poisson(link=identity)) and optim() in R

Following the above question...

And actually, I have an another question now. I could cross validation for glm(poisson, log) to understand Jacobian matrix from log-likelihood by the following program . But I confuse... I can't do with “link=identity”. I wonder how should I change the program for glm(poisson, identity)... Someone, please give me some advice.

x<-c(1,2,3,4)
y<-c(2,3,5,4)
f1<-expression(sum(y-exp(a+b*x)))#simultaneous equations1
f2<-expression(sum(x*(y-exp(a+b*x)))) #simultaneous equations2
g11<-expression(-sum(exp(a+b*x)))#first‐order differentiation for f1 by a
g12<-expression(-sum(x*exp(a+b*x)))#first‐order differentiation for f1 by b
g21<-expression(-sum(x*exp(a+b*x)))#first‐order differentiation for f2 by a
g22<-expression(-sum(x**2*exp(a+b*x)))#first‐order differentiation for f2 by b
a<-1#initial value
b<-0.5#initial value
for(i in 1:10){
m<-matrix(c(a,b),2,1)#a,b to matrix
f<-matrix(c(eval(f1),eval(f2)),2,1)#f1,f2 to matrix
j<-matrix(c(eval(g11), eval(g21), eval(g12), eval(g22)), 2, 2)#Jacobian matrix
m<-m-solve(j)%*%f
a<-m[1]
b<-m[2]
print(sprintf("[%d] (a,b)=(%f,%f)",i,a,b))
}
#cut
#[1] "[7] (a,b)=(0.639265,0.232040)"
#[1] "[8] (a,b)=(0.639265,0.232040)"
#[1] "[9] (a,b)=(0.639265,0.232040)"
#[1] "[10] (a,b)=(0.639265,0.232040)"
(SEs<-sqrt(abs(diag(solve(j)))))
#[1] 0.7317073 0.2445158

#            Estimate Std. Error z value Pr(>|z|)
#(Intercept)   0.6393     0.7317   0.874    0.382
#x             0.2320     0.2445   0.949    0.343


I got the value. But it still didn’t have validation of SE. I don’t know why…

x<-c(1,2,3,4)
y<-c(2,3,5,4)
f<-expression((exp(-a-b)*(a+b)^2/2)*(exp(-a-2*b)*(a+2*b)^3/6)*(exp(-a-3*b)*(a+3*b)^5/120)*(exp(-a-4*b)*(a+4*b)^4/24))
f1<-D(f,"a")
f2<-D(f,"b")
g11<-D(f1,"a")
g12<-D(f1,"b")
g21<-D(f2,"a")
g22<-D(f2,"b")
a<-1.2
b<-0.7
for(i in 1:10){
m<-matrix(c(a,b),2,1)
f<-matrix(c(eval(f1),eval(f2)),2,1)
j<-matrix(c(eval(g11), eval(g21), eval(g12), eval(g22)), 2, 2)
m<-m-solve(j)%*%f
a<-m[1]
b<-m[2]
print(sprintf("[%d] (a,b)=(%f,%f)", i, a, b))
}
#cut
#[1] "[7] (a,b)=(1.278347,0.888661)"
#[1] "[8] (a,b)=(1.278347,0.888661)"
#[1] "[9] (a,b)=(1.278347,0.888661)"
#[1] "[10] (a,b)=(1.278347,0.888661)"
(SEs<-sqrt(abs(diag(solve(j)))))
#[1] 50.43801 20.72464