Multivariate Statistical Modeling Based on Generalized Linear Models

#spline model

#Multivariate Statistical Modeling Based on Generalized Linear Models p.176

library(dplyr)

pis=c(32,104,206,186,102,2,12,28,28,31)/100

x=c(1:length(pis))

h=c()

for(j in 1:(length(x)-1)){

h=c(h,x[j+1]-x[j])

}

D=array(0,dim=c(length(x)-2,length(x)))

C=array(0,dim=c(length(x)-2,length(x)-2))

for(j in 1:nrow(D)){

vec=c(1/h[j],-(1/h[j]+1/h[j+1]),1/h[j+1])

D[j,j:(j+2)]=vec

if(j==1){

C[j,1:2]=c(2*(h[1]+h[2]),h[2])

}else{

if(j==nrow(C)){

C[j,(j-1):j]=c(h[j],2*(h[j]+h[j+1]))

}else{

C[j,(j-1):(j+1)]=c(h[j],2*(h[j]+h[j+1]),h[j+1])

}}

}

K=t(D)%*%solve(C)%*%D

lam=1

f=rep(0.01,length(pis))

sig=rep(1,length(f))

ite=1000;eta=1

for(j in 1:ite){

f_pre=f

D=diag(c(exp(f)))

W=diag(c(diag(D^2)/sig))

y=eta*(pis-exp(f))/diag(D)+f

f=solve(W+lam*K)%*%W%*%y

print(sum(abs(f-f_pre)))

}

S_lam=solve(diag(rep(1,nrow(K)))+lam*K)

GCV=mean(((pis-exp(f))^2/sig)/((1-sum(diag(S_lam))/length(x))^2))

plot(c(1:length(f_pre)),pis,xlim=c(1,length(f_pre)),ylim=c(min(f_pre,f),max(f_pre,f)),col=2,ylab="predict",xlab="index")

par(new=T)

plot(c(1:length(f_pre)),exp(f),xlim=c(1,length(f_pre)),ylim=c(min(f_pre,f),max(f_pre,f)),col=3,ylab="predict",xlab="index")