####################################
#  第4章 モデリング
####################################
この章のRコードは不完全なものがある．後日修正予定
par(mar=c(2,2,1,1)+0.1)


# 多項式回帰
x <- seq(0,1,length=21)
y <- c(0.125, 0.156, 0.193, -0.032, -0.075, -0.064, 0.006, -0.135, 0.105, 0.131, 0.154, 0.114, -0.094, 0.215, 0.035, 0.327, 0.061, 0.383, 0.357, 0.605, 0.499)
#
plot( x,y,pch=16,ylim=c(-0.2,0.6),xlim=c(0,1) )
pol <- polreg(y,7)　# 多項式回帰 0次～6次

#
# このコードは未完成
#
t <- seq(0,21,length=105)
tt <- seq(0,105,length=105)
ey <- rep(0,105)
lag <- 6
s <- 2:lag
z <- pol$coef[[lag]]
for (i in tt){
 ey[i] <- z[1] 
 for (j in s) ey[i] = ey[i] + z[j]*t[i]^(j-1)
}
par(new=T)
plot(t,ey,type="l",lwd=2,ylim=c(-0.2,0.6),xlim=c(0,21),xaxt="n") 


#
# 回帰モデルAIC
#
x <- c(0,1,1,1,2,2,2,3)
y <- c(154.9,119.7,128.3,151.9,119.7,88.9,122.4,90.8)
plot(x,y,pch=16,ylim=c(80,160))

x<- c(-1,0,1,2,3,4,5,6,7,8,9)
y<-c(0.678301,0.006229,0.002587,0.000922,0.000833,0.000737,0.000688,0.000650,0.000622,0.000607,0.000599)
plot(x,y,pch=16,ylim=c(0.0,0.008))


###　確率分布モデルの表示
#　講義で説明していない
#　　TSSSパッケージ
  gdensty(xmin=-4, xmax=4)     # normal distribution
  cdensty(xmin=-4, xmax=4)     # Cauchy distribution
  pdensty(shape=2, xmin=-4, xmax=4)     # Pearson distribution
  exdensty(xmin=0, xmax=8)     # exponential distribution
  exdensty(xmin=-4, xmax=4)    # exponential distribution
  chi2densty(df=3, xmin=0, xmax=8)     # Chi-square distribution
  dexdensty(xmin=-4, xmax=2)   # double exponential distribution
  udensty(xmin=0, xmax=1)     # uniform distribution


###  正規乱数の生成
# Gauss noise: 
z <- simssm(trend=NULL,arcoef=NULL, ar=NULL, seasonal.order=0, seasonal=NULL, 
         tau1=NULL, tau2=NULL, tau3=NULL, n=20, ns=1, init=1992092521, 
         sigma2=1.0, plot=TRUE)
hist(z$y,breaks=seq(-3,3,0.25))


########
stripchart(z$y,pch=1,cex=1.5,xlim=c(-2,2)) 


#############################################
#Computation of Kullback-Leibler Information by numerical integration
# この件は講義で説明していない
# g:gauss, f:gauss
  klinfo(1, c(0, 1), 1, c(0.1, 1.5), 8)

# g:gauss, f:cauchy
  klinfo(1, c(0, 1), 2, c(0, 1), 8)


##############################################
# Box-Cox transformation  # これは第4回の講義用
##############################################
# read subspot data and Whard data
sunspot <- as.ts(read.csv("sunspot_new.csv"))
whard <- as.ts(read.csv("whard_new.csv"))

  data(Sunspot)   # Sun spot number data
  boxcox(Sunspot)

  data(Whard)     # Wholesale hardware data
  boxcox(Whard)

# Box-Cox変換のAICの図示
x <- boxcox(Whard)
plot(x$aic.z,col="red")

# Box-Cox変換の計算結果の表示（AIC,対数尤度，平均，分散
x <- boxcox(Whard)
x