seizure<- read.csv('C:/Users/UX303LB/Desktop/class/R/seizure.csv')
head(seizure)
glm(seizure$y~seizure$trt+seizure$age,offset=seizure$ltime,family=poisson)

#第一題
y<-(seizure$y)
X<-cbind(1,seizure$trt,seizure$age)
head(X)
time<-exp(seizure$ltime)

f3<-function(beta){
  mu<-exp(X%*%beta)*time
  ymu<-y-mu
  gradient<-t(X)%*%ymu
  hessian<--t(X)%*%diag(as.vector(mu))%*%X
  return(list(gradient, hessian))
}

newton <- function(f3, x0, tol = 1e-9, n.max = 100) {
  x <- x0
  f3.x <- f3(x)
  n <- 0
  while ((max(abs(f3.x[[1]])) > tol) & (n < n.max)) {
    x <- x - solve.default(f3.x[[2]], tol= 1e-19)%*% f3.x[[1]]
    f3.x <- f3(x)
    n <- n + 1


  }
  if (n == n.max) {
    cat('newton failed to converge\n')
  } else {
    return(x)
  }
}

newton(f3,c(0, 0, 0),1e-9)

#第二題
result<-glm(seizure$y~seizure$trt+seizure$age,offset=seizure$ltime,family=poisson)
b<-newton(f3,c(0, 0, 0),1e-9)
vcov(result)
solve(-hessian)
solve(-f3(b)[[2]], tol= 1e-19)


#第三題
beta<-newton(f3,c(0, 0, 0),1e-9)
result<-glm(seizure$y~seizure$trt+seizure$age,offset=seizure$ltime,family=poisson)
logLik(result)
loglike<-sum((-exp(X%*%beta)%*%time)+(y%*%X%*%beta)+(y%*%(seizure$ltime))-(log(factorial(y))))

a<--exp(X%*%beta)%*%time
b<-y%*%X%*%beta
c<-y%*%(seizure$ltime)
d<--log(factorial(y))
loglike<-sum(a+diag(b)+diag(c)+d)

loglike<-sum(-mu+diag(y%*%log(mu))-log(factorial(y)))
loglike<-sum(-mu+diag(y%*%log(mu))-log(factorial(y)))