認知情報解析演習 DL-V2 ch4

# p136
embed.forwd <- function(W, idx){
  return(W[idx,])
}

h1 <- embed.forwd(W,idx1)
h2 <- embed.forwd(W,idx2)
h = (h1 + h2)/2

embed.dot.forwd <- function(W, h, idx){
  return(O=rowSums(W[idx,]*h))
}

s = embed.dot.forwd(network$W.out,h,batch.y)

sigmoid.forwd <- function(O){
  return(1/(1+exp(-O)))
}

sigmoidWL.forwd <- function(O,target){
  y = 1/(1+exp(-O)) 
  loss=-sum(target*log(y)+(1-target)*log(1 - y))
  return(loss)
}

sigmoidWL.backwd <- function(O,target,dout=1){
  y = 1/(1+exp(-O)) 
  dW = (y - target)*dout
  return(dW)
}

embed.dot.backwd <- function(W,h,idx,dout){
  dh = dout * W[idx,]
  dW = dout * h
  return(list(dh = dh, dW = dW))
}

embed.backwd <- function(W, idx, dout){
  dW = matrix(0,nrow(W),ncol(W))
  for (i.idx in 1:length(idx)){
    dW[idx[i.idx], ] = dW[idx[i.idx], ] + dout[i.idx, ]
  }
  return(dW)
}

# 158

txt2corpus <- function(txt){
  txt = tolower(txt)
  txt = gsub('[.]', ' .',txt)
  words = unlist(strsplit(txt, " "))
  uniq.words = unique(words)
  n.uniq = length(uniq.words)
  n.words = length(words)
  corpus = rep(0,n.words)
  corpus = match(words,uniq.words)
  return(corpus)
}

corp2contxt <- function(corpus){
  len.corp = length(corpus)
  target = corpus[2:(len.corp-1)]
  col1 = corpus[1:(len.corp-2)]
  col2 = corpus[3:len.corp]
  context = cbind(col1,col2)
  return(list(context=context,target = target))
}

init.W2V <- function(n.uniq,size.hidden){
  W.in = matrix(rnorm(n.uniq*size.hidden),nrow=n.uniq)*0.01
  W.out = matrix(rnorm(n.uniq*size.hidden),nrow=n.uniq)*0.01
  return(list(W.in = W.in, W.out= W.out))
}


### running example
network<-init.W2V(8,5)
n.batch = 3;
txt = "You say goodbye and I say hello.."
corp = txt2corpus(txt)
corp = c(8,8,corp)
dat<-corp2contxt(corp)
n.data =  length(dat$target)

n.iter = 3000;
lambda=0.1;
loss = rep(0, n.iter)
for (i.iter in 1:n.iter){
  samp.idx = sample(1:n.data,n.batch)
  batch.c1 = dat$context[samp.idx,1]
  batch.c2 = dat$context[samp.idx,2]
  batch.y = dat$target[samp.idx]
  h1 <- embed.forwd(network$W.in,batch.c1)
  h2 <- embed.forwd(network$W.in,batch.c2)
  h = 0.5 * (h1 + h2)
  # positive only
  s = embed.dot.forwd(network$W.out,h,batch.y)
  z = sigmoidWL.forwd(s,1)
  loss[i.iter] = z
  ds = sigmoidWL.backwd(s, 1, 1)
  dE = embed.dot.backwd(network$W.out,h, batch.y, ds)
  dh = dE$dh*0.5
  dIn1 = embed.backwd(network$W.in,dat$context[samp.idx,1], dh)
  dIn2 = embed.backwd(network$W.in,dat$context[samp.idx,2], dh)
  network$W.out[batch.y,] = network$W.out[batch.y,] - lambda*dE$dW
  network$W.in = network$W.in - lambda*dIn1
  network$W.in = network$W.in - lambda*dIn2
}
plot(loss, type='l')

samp.idx = c(2:6,8,9,1)
batch.c1 = dat$context[samp.idx,1]
batch.c2 = dat$context[samp.idx,2]
batch.y = dat$target[samp.idx]
h1 <- embed.forwd(network$W.in,batch.c1)
h2 <- embed.forwd(network$W.in,batch.c2)
h = 0.5 * (h1 + h2)
s = MatMult.forwd(h,t(network$W.out))
z = sigmoid.forwd(s)
res=cbind(z,batch.y)
par(mfrow=c(8,1))
for (i in 1:8){
  col.spec = rep("black",8)
  col.spec[i]="orange"
  barplot(res[i, 1:8],col=col.spec)
}
# with negative sampling
# p158

unigram.sampler <- function(corpus, target, power, sample.size){
  neg.corpus <- corpus[-which(corpus == target)]
  uniq.word = unique(neg.corpus)
  n.word = length(neg.corpus)
  prob = (as.vector(table(neg.corpus))/n.word)^power
  prob = prob/sum(prob)
  sample.idx = sample(uniq.word, sample.size, prob = prob)
  return(sort(sample.idx))
}

get.neg.sample <- function(corpus, target, power, sample.size){
  sample.id = matrix(0, nrow = length(target), ncol = sample.size)
  for (i.targ in 1:length(target)){
    sample.id[i.targ,] = unigram.sampler(corpus, target[i.targ], power, sample.size)
  }
  return(sample.id)
}

init.W2V <- function(n.uniq,size.hidden){
  W.in = matrix(rnorm(n.uniq*size.hidden),nrow=n.uniq)*0.01
  W.out = matrix(rnorm(n.uniq*size.hidden),nrow=n.uniq)*0.01
  return(list(W.in = W.in, W.out= W.out))
}

network<-init.W2V(8,5)
n.batch = 3;
txt = "You say goodbye and I say hello.."
corp = txt2corpus(txt)
corp = c(8,8,corp)
dat<-corp2contxt(corp)
n.data =  length(dat$target)

n.iter = 3000;
lambda=0.1;
loss = rep(0, n.iter)
sample.size = 2
for (i.iter in 1:n.iter){
  samp.idx = sample(1:n.data,n.batch)
  batch.c1 = dat$context[samp.idx,1]
  batch.c2 = dat$context[samp.idx,2]
  batch.y = dat$target[samp.idx]
  h1 <- embed.forwd(network$W.in,batch.c1)
  h2 <- embed.forwd(network$W.in,batch.c2)
  h = 0.5 * (h1 + h2)
  # positive
  s = embed.dot.forwd(network$W.out,h,batch.y)
  z = sigmoidWL.forwd(s,1)
  loss[i.iter] = z
  ds = sigmoidWL.backwd(s, 1, 1)
  dE = embed.dot.backwd(network$W.out,h, batch.y, ds)
  dh = dE$dh*0.5
  dIn1 = embed.backwd(network$W.in,dat$context[samp.idx,1], dh)
  dIn2 = embed.backwd(network$W.in,dat$context[samp.idx,2], dh)
  network$W.out[batch.y,] = network$W.out[batch.y,] - lambda*dE$dW
  network$W.in = network$W.in - lambda*dIn1
  network$W.in = network$W.in - lambda*dIn2
  # negative
  neg.set <- get.neg.sample(corp,batch.y, 0.75, sample.size)
  for (i.neg in 1:sample.size){
      s = embed.dot.forwd(network$W.out,h,neg.set[,i.neg])
      z = sigmoidWL.forwd(s,0)
      loss[i.iter] = loss[i.iter] + z
      ds = sigmoidWL.backwd(s, 0, dout=1)
      dE = embed.dot.backwd(network$W.out,h,neg.set[,i.neg], ds)
      dh = dE$dh*0.5
      dIn1 = embed.backwd(network$W.in, dat$context[samp.idx,1],dh)
      dIn2 = embed.backwd(network$W.in, dat$context[samp.idx,2],dh)
      network$W.out[neg.set[,i.neg],] = network$W.out[neg.set[,i.neg],] - lambda/sample.size*dE$dW
      network$W.in = network$W.in - lambda*dIn1
      network$W.in = network$W.in - lambda*dIn2
    }
}
plot(loss, type='l')

samp.idx = c(2:6,8,9,1)
batch.c1 = dat$context[samp.idx,1]
batch.c2 = dat$context[samp.idx,2]
batch.y = dat$target[samp.idx]
h1 <- embed.forwd(network$W.in,batch.c1)
h2 <- embed.forwd(network$W.in,batch.c2)
h = 0.5 * (h1 + h2)
s = MatMult.forwd(h,t(network$W.out))
z = sigmoid.forwd(s)
res=cbind(z,batch.y)
par(mfrow=c(8,1))
for (i in 1:8){
  col.spec = rep("black",8)
  col.spec[i]="orange"
  barplot(res[i, 1:8],col=col.spec)
}

par(mfrow=c(4,1))
for (i in 1:4){
  col.spec = rep("black",8)
  col.spec[i]="orange"
  barplot(res[i, 1:8],col=col.spec)
}

院:認知情報解析学演習 ch23 & 24

dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/OrdinalProbitData-1grp-1.csv" )

model = "
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dcat( pr[i,1:nYlevels] )
    pr[i,1] <- pnorm( thresh[1] , mu , 1/sigma^2 )
    for ( k in 2:(nYlevels-1) ) {
      pr[i,k] <- max( 0 ,  pnorm( thresh[ k ] , mu , 1/sigma^2 )- pnorm( thresh[k-1] , mu , 1/sigma^2 ) )
    }
    pr[i,nYlevels] <- 1 - pnorm( thresh[nYlevels-1] , mu , 1/sigma^2 )
  }
  mu ~ dnorm( (1+nYlevels)/2 , 1/(nYlevels)^2 )
  sigma ~ dunif( nYlevels/1000 , nYlevels*10 )
  for ( k in 2:(nYlevels-2) ) {  # 1 and nYlevels-1 are fixed, not stochastic
    thresh[k] ~ dnorm( k+0.5 , 1/2^2 )
  }
}
"
writeLines( model , "model.txt" )

y = dat$Y
Ntotal = length(y)
nYlevels = max(y)
thresh = rep(NA,nYlevels-1)
thresh[1] = 1 + 0.5
thresh[nYlevels-1] = nYlevels-1 + 0.5
dataList = list(y = y , nYlevels = nYlevels,thresh = thresh, Ntotal = Ntotal)

parameters = c( "pr" ,  "thresh","mu","sigma")
jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)

meanT = rowMeans(cbind(mcmcMat[,"thresh[1]"],mcmcMat[,"thresh[2]"],mcmcMat[,"thresh[3]"],mcmcMat[,"thresh[4]"],mcmcMat[,"thresh[5]"],
                   mcmcMat[,"thresh[6]"]))
plot(mcmcMat[,"thresh[1]"],meanT,xlim=c(0.5,7))
points(mcmcMat[,"thresh[2]"],meanT,xlim=c(0.5,7))
points(mcmcMat[,"thresh[3]"],meanT,xlim=c(0.5,7))
points(mcmcMat[,"thresh[4]"],meanT,xlim=c(0.5,7))
points(mcmcMat[,"thresh[5]"],meanT,xlim=c(0.5,7))
points(mcmcMat[,"thresh[6]"],meanT,xlim=c(0.5,7))
HDI.plot(mcmcMat[,"mu"],xlab="mu")
HDI.plot(mcmcMat[,"sigma"],xlab="sigma")


model2 = "
  model {
for ( i in 1:Ntotal ) {
y[i] ~ dcat( pr[i,1:nYlevels] )
pr[i,1] <- pnorm( thresh[1] , mu[x[i]] , 1/sigma[x[i]]^2 )
for ( k in 2:(nYlevels-1) ) {
pr[i,k] <- max( 0 ,  pnorm( thresh[ k ] , mu[x[i]] , 1/sigma[x[i]]^2 )
- pnorm( thresh[k-1] , mu[x[i]] , 1/sigma[x[i]]^2 ) )
}
pr[i,nYlevels] <- 1 - pnorm( thresh[nYlevels-1] , mu[x[i]] , 1/sigma[x[i]]^2 )
}
for ( j in 1:2 ) { # 2 groups
mu[j] ~ dnorm( (1+nYlevels)/2 , 1/(nYlevels)^2 )
sigma[j] ~ dunif( nYlevels/1000 , nYlevels*10 )
}
for ( k in 2:(nYlevels-2) ) {  # 1 and nYlevels-1 are fixed, not stochastic
thresh[k] ~ dnorm( k+0.5 , 1/2^2 )
}
}
"
writeLines( model2 , "model2.txt" )

dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/OrdinalProbitData1.csv")
y = dat$Y
x = as.numeric(dat$X)
Ntotal = length(y)
nYlevels = max(y)
thresh = rep(NA,nYlevels-1)
thresh[1] = 1 + 0.5
thresh[nYlevels-1] = nYlevels-1 + 0.5
dataList = list(y = y , x= x ,nYlevels = nYlevels,thresh = thresh, Ntotal = Ntotal)

parameters = c( "pr" ,  "thresh","mu","sigma")
jagsModel = jags.model( "model2.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)
HDI.plot(mcmcMat[,"mu[1]"])
HDI.plot(mcmcMat[,"mu[2]"])
HDI.plot(mcmcMat[,"mu[1]"]-mcmcMat[,"mu[2]"])


model3 = "
data {
xm  <- mean(x)
xsd <-   sd(x)
for ( i in 1:Ntotal ) {
zx[i] <- ( x[i] - xm ) / xsd
}
}
model {
for ( i in 1:Ntotal ) {
y[i] ~ dcat( pr[i,1:nYlevels] )
pr[i,1] <- pnorm( thresh[1] , mu[i] , 1/sigma^2 )
for ( k in 2:(nYlevels-1) ) {
pr[i,k] <- max( 0 ,  pnorm( thresh[ k ] , mu[i] , 1/sigma^2 )
- pnorm( thresh[k-1] , mu[i] , 1/sigma^2 ) )
}
pr[i,nYlevels] <- 1 - pnorm( thresh[nYlevels-1] , mu[i] , 1/sigma^2 )
mu[i] <- zbeta0 + sum( zbeta * zx[i] )
}
# Priors vague on standardized scale:
zbeta0 ~ dnorm( (1+nYlevels)/2 , 1/(nYlevels)^2 )
zbeta ~ dnorm( 0 , 1/(nYlevels)^2 )

zsigma ~ dunif( nYlevels/1000 , nYlevels*10 )
# Transform to original scale:
beta <- ( zbeta / xsd )
beta0 <- zbeta0  - sum( zbeta * xm / xsd )
sigma <- zsigma
for ( k in 2:(nYlevels-2) ) {  # 1 and nYlevels-1 are fixed
thresh[k] ~ dnorm( k+0.5 , 1/2^2 )
}
}
"
writeLines( model3, con="model3.txt" )

dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/OrdinalProbitData-LinReg-2.csv")
y = dat$Y
x = as.numeric(dat$X)
Ntotal = length(y)
Nx =1
nYlevels = max(y)
thresh = rep(NA,nYlevels-1)
thresh[1] = 1 + 0.5
thresh[nYlevels-1] = nYlevels-1 + 0.5
dataList = list(y = y , Nx= 1,x= x ,nYlevels = nYlevels,thresh = thresh, Ntotal = Ntotal)

parameters = c("thresh","mu","sigma","beta","beta0")
jagsModel = jags.model( "model3.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)
plot(dat$X,dat$Y,pch=20,cex=3)
for (i.plot in 1:100){
  abline(a=mcmcMat[i.plot,"beta0"],b=mcmcMat[i.plot,"beta"],col="orange",lwd=2)
}
HDI.plot(mcmcMat[,"beta"])
#plot(1,1,type='n',xlim=c(0,6),ylim=c(0,6))
#abline(a=0,b=1,lwd=2)
#x.temp = seq(0,6,length.out=100)
#y = dnorm(x.temp,mean=3,sd=1)
#lines(x=-y*3+1,y=x.temp)
#lines(x=c(-1,3),y=c(3,3))
#lines(x=c(-1,0.5),y=c(0.5,0.5),lty=2,col="red")
#lines(x=c(-1,1.5),y=c(1.5,1.5),lty=2,col="red")
#lines(x=c(-1,2.5),y=c(2.5,2.5),lty=2,col="red")
#lines(x=c(-1,3.5),y=c(3.5,3.5),lty=2,col="red")
#lines(x=c(-1,4.5),y=c(4.5,4.5),lty=2,col="red")
#lines(x=c(-1,5.5),y=c(5.5,5.5),lty=2,col="red")

model4="
model {
  for ( i in 1:Ncell ) {
    y[i] ~ dpois( lambda[i] )
    lambda[i] <- exp( a0 + a1[x1[i]] + a2[x2[i]] + a1a2[x1[i],x2[i]] )
  }
  a0 ~ dnorm( yLogMean , 1/(yLogSD*2)^2 )
  for ( j1 in 1:Nx1Lvl ) { a1[j1] ~ dnorm( 0.0 , 1/a1SD^2 ) }
  a1SD ~ dgamma(agammaShRa[1],agammaShRa[2])
  for ( j2 in 1:Nx2Lvl ) { a2[j2] ~ dnorm( 0.0 , 1/a2SD^2 ) }
  a2SD ~ dgamma(agammaShRa[1],agammaShRa[2])
  for ( j1 in 1:Nx1Lvl ) { for ( j2 in 1:Nx2Lvl ) {
    a1a2[j1,j2] ~ dnorm( 0.0 , 1/a1a2SD^2 )
  } }
  a1a2SD ~ dgamma(agammaShRa[1],agammaShRa[2])
  # Convert a0,a1[],a2[],a1a2[,] to sum-to-zero b0,b1[],b2[],b1b2[,] :
  for ( j1 in 1:Nx1Lvl ) { for ( j2 in 1:Nx2Lvl ) {
    m[j1,j2] <- a0 + a1[j1] + a2[j2] + a1a2[j1,j2] # cell means
  } }
  b0 <- mean( m[1:Nx1Lvl,1:Nx2Lvl] )
  for ( j1 in 1:Nx1Lvl ) { b1[j1] <- mean( m[j1,1:Nx2Lvl] ) - b0 }
  for ( j2 in 1:Nx2Lvl ) { b2[j2] <- mean( m[1:Nx1Lvl,j2] ) - b0 }
  for ( j1 in 1:Nx1Lvl ) { for ( j2 in 1:Nx2Lvl ) {
    b1b2[j1,j2] <- m[j1,j2] - ( b0 + b1[j1] + b2[j2] )
  } }
  # Compute predicted proportions:
  for ( j1 in 1:Nx1Lvl ) { for ( j2 in 1:Nx2Lvl ) {
    expm[j1,j2] <- exp(m[j1,j2])
    ppx1x2p[j1,j2] <- expm[j1,j2]/sum(expm[1:Nx1Lvl,1:Nx2Lvl])
  } }
  for ( j1 in 1:Nx1Lvl ) { ppx1p[j1] <- sum(ppx1x2p[j1,1:Nx2Lvl]) }
  for ( j2 in 1:Nx2Lvl ) { ppx2p[j2] <- sum(ppx1x2p[1:Nx1Lvl,j2]) }
}"
writeLines( model4, con="model4.txt" )

dat = read.csv( file="http://peach.l.chiba-u.ac.jp/course_folder/HairEyeColor.csv" )

y=dat$Count
x1=dat$Eye
x2=dat$Hair

x1 = as.numeric(as.factor(x1))
x1levels = levels(as.factor(dat$Eye))
x2 = as.numeric(as.factor(x2))
x2levels =  levels(as.factor(dat$Eye))
Nx1Lvl = length(unique(x1))
Nx2Lvl = length(unique(x2))
Ncell = length(y) # number of rows of data, not sum(y)

yLogMean = log(sum(y)/(Nx1Lvl*Nx2Lvl))
yLogSD = log(sd(c(rep(0,Ncell-1),sum(y))))
MeanSD2gamma <- function( mean, sd ) {
  shape = mean^2 / sd^2
  rate = mean / sd^2
  return(data.frame(shape,rate))
}

ModeSD2gamma <- function( mode, sd ) {
  rate = ( mode + sqrt( mode^2 + 4 * sd^2 ) )/( 2 * sd^2 )
  shape = 1 + mode * rate
  return(data.frame(shape,rate))
}
temp=ModeSD2gamma(mode=yLogSD , sd=2*yLogSD)
agammaShRa = unlist(temp )

data_list = list(
  y = y ,x1 = x1 , x2 = x2 ,
  Ncell = Ncell , Nx1Lvl = Nx1Lvl , Nx2Lvl = Nx2Lvl ,
  yLogMean = yLogMean , yLogSD = yLogSD ,agammaShRa = agammaShRa
)

jagsModel =jags.model("model4.txt", data = data_list, n.chains = 3, n.adapt = 500)
update(jagsModel, 500)
codaSamples = coda.samples(jagsModel, variable.names = c("b0", "b1", "b2", "b1b2", "ppx1p", "ppx2p", "ppx1x2p"),
                           n.iter = ((10000*1)/1), n.adapt = 500)
mcmcMat<-as.matrix(codaSamples)
EyeBlueHairBlack <- mcmcMat[,"ppx1x2p[1,1]"]
HDI.plot(EyeBlueHairBlack)


EyeBlue_vs_EyeBrown_at_HairBlack <- mcmcMat[,"b1b2[1,1]"] - mcmcMat[,"b1b2[2,1]"] - mcmcMat[,"b1[2]"] + mcmcMat[,"b1[1]"]
HDI.plot(EyeBlue_vs_EyeBrown_at_HairBlack)

EyeBlue_vs_EyeBrown_at_Blond <- mcmcMat[,"b1b2[1,2]"] - mcmcMat[,"b1b2[2,2]"] - mcmcMat[,"b1[2]"] + mcmcMat[,"b1[2]"]
HDI.plot(EyeBlue_vs_EyeBrown_at_Blond)

diff <- EyeBlue_vs_EyeBrown_at_HairBlack - EyeBlue_vs_EyeBrown_at_Blond
HDI.plot(diff)

基礎自習B02

dec2bin<-function(num, digits=8) {
  bin=c()
  if (num==0){
    bin=0
  } else {
    while(num!=0){
      rem=num%%2
     num=num%/%2
      bin=c(rem,bin)
    }
  }
  if (length(bin) < digits){
    res=matrix(0,nrow=1,ncol=digits)
    res[(digits-length(bin)+1):digits]=bin
  } else {res=bin}
  return(res)
}

データ解析基礎論B Loglineak Analysis

freq<-c(33,37,7,23)
pref<-factor(c('soba','udon','soba','udon'))
region<-factor(c('east','east','west','west'))
dat<-data.frame(pref,region,freq)
dat.table=table(pref,region)
dat.table[cbind(pref,region)]<-freq

dat.loglinCE_A<-loglin(dat.table, list(1), fit=T,param=T)
dat.loglinCE_B<-loglin(dat.table,list(2), fit=T,param=T)
dat.loglinIND<-loglin(dat.table,list(1,2), fit=T,param=T)
dat.loglinSAT<-loglin(dat.table,list(c(1,2)), fit=T,param=T)


freq<-c(9,5,2,4,16,10,26,28)
gener<-factor(c(rep('female',4),c(rep('male',4))))
affil<-factor(rep(c('L','L','E','E'),2))
circle<-factor(rep(c('tennis','astro'),4))
dat<-data.frame(gener,affil,circle,freq)

dat.table<-table(gender,affil,circle)
dat.table[cbind(gender,affil,circle)]<-freq

dat.loglin2<-loglin(dat.table,list(1), fit=T,param=T)
dat.loglin3<-loglin(dat.table,list(1,3), fit=T,param=T)
dat.loglin4<-loglin(dat.table,list(1,2,3), fit=T,param=T)
dat.loglin5<-loglin(dat.table,list(c(1,3),2), fit=T,param=T)
dat.loglin6<-loglin(dat.table,list(c(1,3),c(1,2)), fit=T,param=T)
dat.loglin7<-loglin(dat.table,list(c(1,3),c(1,2),c(2,3)), fit=T,param=T)
dat.loglin8<-loglin(dat.table,list(c(1,2,3)), fit=T,param=T)

source('http://peach.l.chiba-u.ac.jp/course_folder/cuUtil02.R')