Disclaimer

このcourselogにあるコードは、主に学部生・博士課程前期向けの講義・演習・実習で例として提示しているもので、原則直感的に分かりやすいように書いているつもりです。例によってはとても非効率なものもあるでしょうし、「やっつけ」で書いているものあります。また、普段はMATLABを使用していますので、変な癖がでているかもしれません。これらの例を使用・参考にする場合はそれを踏まえてたうえで使用・参考にして下さい。
卒業論文に関する資料:[2015PDF] [word template] [latex template] [表紙] [レポートの書き方] [引用文献など]
A Brief Guide to R (Rのコマンドの手引きおよび例)はこちらをどうぞ

データ解析基礎論B W03 R graphics

library(tidyverse)


# CLT
NperSample = 10
SampleSize = 300000

runif(NperSample * SampleSize) %>%
  matrix(nrow=NperSample) %>%
     colMeans() %>% tibble(sample.mean = .) -> means

ggplot(means,aes(x = sample.mean, y = ..density..)) +
  geom_histogram(bins=200) +
    geom_density(colour = "orange",size=2)

ggplot(means,aes(x = sample.mean, y = ..density..)) +
  geom_histogram(bins=200) +
  geom_line(stat = "density", colour = "orange",size=2)

runif(NperSample * SampleSize) %>%
  matrix(nrow=NperSample) %>%
  colMeans() %>% tibble(sample.mean = .) %>%
  ggplot(., aes(x = sample.mean, y = ..density..)) +
    geom_histogram(bins=100,colour = "grey20") +
    geom_line(stat = "density", colour = "skyblue",size=2)


dat <- read.csv("http://www.matsuka.info/data_folder/sampleData2013.txt")
dt <- as_tibble(dat)




ggplot(dt, aes(x = Hworked, y = nbooks)) +
  geom_point(size = 3)

ggplot(dt) +
  geom_point(aes(x = Hworked, y = nbooks, color = grade),size = 3)

ggplot(dt) +
  geom_point(aes(x = Hworked, y = nbooks, shape = grade),size = 5)

ggplot(dt) +
  geom_point(aes(x = Hworked, y = nbooks),size = 5) +
  facet_wrap(~ grade, nrow = 1)

ggplot(dt) +
  geom_smooth(aes(x = Hworked, y = nbooks))

ggplot(dt) +
  geom_smooth(aes(x = Hworked, y = nbooks, linetype = grade))

ggplot(dt) +
  geom_smooth(aes(x = Hworked, y = nbooks)) +
    facet_wrap(~ grade, nrow = 4)

ggplot(dt) +
  geom_smooth(aes(x = Hworked, y = nbooks)) +
  geom_point(aes(x = Hworked, y = nbooks), size = 4)

ggplot(dt) +
  geom_smooth(aes(x = Hworked, y = nbooks), colour = "black", se = FALSE) +
  geom_point(aes(x = Hworked, y = nbooks, color = grade), size = 4)

ggplot(dt) +
  geom_smooth(aes(x = Hworked, y = nbooks, color = grade), se = FALSE) +
  geom_point(aes(x = Hworked, y = nbooks, color = grade), size = 4)


plot1 <- ggplot(dt) +
  geom_smooth(aes(x = Hworked, y = nbooks, color = grade), se = FALSE) +
  geom_point(aes(x = Hworked, y = nbooks, color = grade), size = 4)
plot1 + xlab("Hours worked") + ylab("Number of books read")

plot1 + xlab("Hours worked") +  ylab("Number of books read") +
  theme(axis.title.x = element_text(face = "italic",size = 14, colour = "navy"),
        axis.title.y = element_text(face = "bold",size = 10, colour = "darkgreen"))

ggplot(filter(dt, affil == "LA")) +
  geom_point(aes(x = Hworked, y = nbooks, color = grade), size = 4)


dt$grade <- fct_relevel(dt$grade, "FR","SP","JR","SR")
group_by(dt, grade) %>% summarize(ave.books = mean(nbooks, na.rm = T)) %>%
  ggplot() + geom_bar(aes(x = grade, y = ave.books), stat = "identity")

group_by(dt, grade) %>% summarize(ave.books = mean(nbooks, na.rm = T)) %>%
  ggplot() + geom_bar(aes(x = grade, y = ave.books), stat = "identity")

group_by(dt, grade) %>% summarize(ave.books = mean(nbooks, na.rm = T),
                                  se = sd(nbooks, na.rm =T)/n()) %>%
ggplot(aes(x = grade, y = ave.books)) +
  geom_bar(stat = "identity", fill = "grey70") +
  geom_errorbar(aes(ymin = ave.books - se, ymax = ave.books +se), width = 0.2) +
  ylab("Average # books read")

ggplot(dt,aes(x = Hworked, y = nbooks)) +
  stat_density2d(aes(colour =..level..)) +
  geom_point()

ggplot(dt,aes(x = Hworked, y = nbooks)) +
  stat_density2d(aes(alpha =..density..), geom="tile",contour=F) +
 geom_point(alpha =0.4)


ggplot(dt) +
  stat_summary(aes(x = grade, y = nbooks),
               fun.y = mean,
               fun.ymin = function(x) mean(x) - sd(x),
               fun.ymax = function(x) mean(x) + sd(x))

ggplot(dt) +
  geom_boxplot(aes(x = grade, y = nbooks))
ggplot(dt) +
  geom_boxplot(aes(x = grade, y = nbooks)) +
  coord_flip()

dat <- read.csv("http://www.matsuka.info/data_folder/datWA01.txt")
dt <- as_tibble(dat)
dt.lm <- lm(h~shoesize, dt)
cfs <- coef(dt.lm)
ggplot(dt, aes(x = shoesize, y = h)) +
  geom_point() +
  geom_abline(intercept = cfs[1], slope = cfs[2], col = "red") +
  geom_text( x= 22, y =175, aes(label = paste("r^2  =",round(summary(dt.lm)$r.squared,3))))

広域システム特別講義II 10.11の講義

10.11ABの講義は家族が急病の為、休講とさせてください。直前の連絡で申し訳ありません。
次回の講義・実習は10.25日です。
補講日は次回、皆さんの都合で調整しましょう。
補講日の候補は11.01、12.27、 および、東大が指定している補講日である、12.25(午後)、01.17(午前)、01.20、01.21です。

Posted in UT

データ解析基礎論B W02

install.packages("tidyverse")
library(tidyverse)

random.number <- rnorm(1000)
mean(random.number)
mean(random.number <- rnorm(1000))

rnorm(1000) %>% mean()

# CLT
NperSample = 10
SampleSize = 300000

# "traditional"
random.number <- runif(NperSample * SampleSize)
dat <- matrix(random.number, nrow=NperSample)
means <- colMeans(dat)
dens <- density(means)
hist(means, breaks = 100, probability = T, main = "Distributionf of Means")
lines(dens, lwd = 3, col = "orange")


runif(NperSample * SampleSize) %>% 
  matrix(nrow=NperSample) %>%
     colMeans() -> means
hist(means, breaks=100,probability = T, main = "Distributionf of Means")
means %>% density() %>% lines(,lwd=3,col='orange')


histWdensity <- function(dat, nbreaks=30, main.txt){
  dens <- density(dat)
  hist(dat, breaks = nbreaks, probability = T, main = main.txt)
  lines(dens, lwd = 3, col = "orange")
}

runif(NperSample * SampleSize) %>% 
  matrix(nrow=NperSample) %>%
  colMeans() %>% 
    histWdensity(nbreaks = 100, main.txt = "Distributionf of Means")

dat <- read.csv("http://www.matsuka.info/data_folder/sampleData2013.txt")
dt <- as_tibble(dat)
dt.la <- filter(dt, affil == "LA")

dt.la2 <- filter(dt, affil == "LA" & grade == "SP")
dt.laNot2 <- filter(dt, affil == "LA" & grade != "SP")

dt.GB <- select(dt, grade, nbooks)
dtWOgender <- select(dt, -gender)

dt.arranged <- arrange(dt, affil, grade)

dt.weekly <- mutate(dt,nbooksWeek = nbooks / 52)

dt.atLeastOneBook <- mutate(dt, atleastOneBook = (nbooks/52) >= 1)
dt.atLeastOneBook <- mutate(dt, atleastOneBook = (nbooks/12) >= 1)

dt.BWindex = mutate(dt, nbooksWeek = nbooks / 52, 
                        idx = nbooksWeek / (12*7-Hworked))

dt.byGrade <- group_by(dt, grade)
summarize(dt.byGrade, ave.books = mean(nbooks,na.rm = TRUE), 
                      ave.Hworked = mean(Hworked, na.rm = TRUE))

dt.byGrAf <- group_by(dt, grade, affil)
dt.summ <- summarize(dt.byGrAf, ave.books = mean(nbooks,na.rm = TRUE),
                     ave.Hworked = mean(Hworked, na.rm = TRUE), N = n())

dt.summ2 <- dt %>% 
  group_by(grade, affil) %>% 
  summarize(ave.books = mean(nbooks,na.rm = TRUE), 
            ave.Hworked = mean(Hworked, na.rm = TRUE), 
            N = n()) %>% filter(N > 2) %>% arrange(desc(ave.books))

plot(x = dt.summ2$ave.books, y = dt.summ2$ave.Hworked, pch=20, cex = 3,
  xlab = "Ave. # books read",ylab = "Ave hours worked")

dt <- read_csv("http://www.matsuka.info/data_folder/sampleData2013.txt",
  col_names = TRUE)

dt.summ3 <- dt %>% 
  group_by(grade, gender) %>% 
  summarize(ave.books = mean(nbooks,na.rm = TRUE), 
    ave.Hworked = mean(Hworked, na.rm = TRUE)) 

dt.summ3G <- dt.summ3 %>% gather('ave.books', 'ave.Hworked', 
  key = 'ave.values', value = "BorW")

dt.summ3SformG <- spread(dt.summ3G, key = ave.values, value =BorW)

dt.sumLA <- dt %>% filter(affil=="LA") %>% group_by(grade) %>% 
   summarize(ave.books = mean(nbooks))

toeic <- tibble(
  grade = factor(c("SP", "JR")),
  score = c(800,830),
)

new.dt1 <- dt.sumLA %>% inner_join(toeic, by = "grade")

dt.sumLA <- add_row(dt.sumLA, grade = c("MS"), ave.books = (13))
toeic2 <- tibble(
  grade = factor(c("SP", "JR","DR")),
  score = c(800,830,900),
)
new.dt3 <- full_join(dt.sumLA, toeic2)

new.dt4 <- left_join(dt.sumLA, toeic2)
new.dt5 <- right_join(dt.sumLA, toeic2)

広域システム特別講義II 課題1

課題1.1 勾配法を用いた最小化

z= \frac{1}{20}x^2+y^2
を最小化するxとyを慣性を含む勾配法を用いて求めてください。

fun01 = function(x,y){
  z = 1/20*x^2 + y^2
  return(z)
}


結果の可視化は以下のようなものでcontourを作図して:

install.packages("plot3D")
library(plot3D)
x = seq(-10,10,0.2)
y = seq(-5,5,0.2)
M = mesh(x,y)
Z = as.vector(1/20*M$x^2)+as.vector(M$y^2) 
Z.mesh = matrix(Z,nrow(M$x))
contour(x, y, Z.mesh, drawlabels = F, nlevels=40)

その後linesで更新履歴を追加してみると良いかもしれません。

# gdは勾配法でのx、yの更新履歴
lines(gd, type='o', col = 'green', pch=20)

# gdmは慣性ありの勾配法でのx、yの更新履歴
lines(gdm, type='o', col = 'blue', pch=20)

課題1.2 multi-armed bandit problem
wikipediaの記事
何種類かバージョンはありますが、この課題では各スロットマシーンの真のリワードは、mean = 0、sd = 1の正規分布に従うこととしてください。ただし、ある特定の回で実際に得られるリワードは真のリワード、プラス正規分布にしたがうノイズがあるとします(ここでも、mean = 0、sd = 1としましょう)。
真のリワードは不明で、1000回スロットマシーンを引いた時の総リワードを最大化することを目的とします。
1つ目の方策はgreedy法で、リワードの推定値が最も高いスロットマシーンを選択するといったものです。
2つ目の方策は基本的にはgreedy法で、ある確率epsilonでランダムにスロットマシーンを選択するものです。
各方策につき複数回検証しましょう(M回)。

# スロットマシーンの数を5としました(実際には幾つでも結構です)
N = 5

# スロットマシーンを引く回数を1000としました(実際には幾つでも結構です)
nTrial = 1000

# 真のリワードを生成します。
reward.true = rnorm(N, mean=0, sd=1)
> reward.true
[1] -0.2822860  0.5645874 -0.1968128  0.5430834 -0.3696859

# リワードの推定値を初期化します。
reward.est = rep(0, N)

# リワードの累積和を初期化します。
reward.cum = rep(0, N)

# 各スロットマシーンを引いた回数を初期化します。
sm.count = rep(0, N)

# リワードの履歴を初期化します。
reward.hist = rep(0, nTrial)

# greedy法で、どのスロットマシーンを引くか選択します。
# reward.estの最大値が複数個ある場合は、それらから1つをランダムで選択します。
max.est = which(max(reward.est) == reward.est)
if (length(max.est) > 1){
  selected = sample(max.est, 1)
} else {selected = max.est}

# 今回は5が選択されました。
> selected
[1] 5

# スロットマシーン5を引いてみます。
# 真のリワードは
# > reward.true[selected]
# [1] -0.3696859
# ですが、今回実際に得られるのはこれにノイズが乗ります。
reward = reward.true[selected] + rnorm(1, mean = 0, sd =1)
> reward
[1] -1.61256

reward.hist[1] = reward
# 繰り返す場合は、reward.hist[i.trial] = reward など 


# リワードの推定値をアップデートします
reward.cum[selected] = reward.cum[selected] + reward
sm.count[selected] = sm.count[selected] + 1
reward.est[selected] = reward.cum[selected] / sm.count[selected]
> reward.est
[1]  0.00000  0.00000  0.00000  0.00000 -1.61256

# 2回目
max.est = which(max(reward.est) == reward.est)
if (length(max.est) > 1){
   selected = sample(max.est, 1)
 } else {selected = max.est}
> selected 
[1] 2

reward = reward.true[selected] + rnorm(1, mean = 0, sd =1)
> reward
[1] 1.497099

reward.cum[selected] = reward.cum[selected] + reward
sm.count[selected] = sm.count[selected] + 1
reward.est[selected] = reward.cum[selected] / sm.count[selected]

> reward.est
[1]  0.000000  1.497099  0.000000  0.000000 -1.612560

# これをnTrial分繰り返します。

# 2つの方策の良さの検証は、特定のreward.trueの値に依存するべきではないので、
# reward.trueの値を変えてみます。これをM回繰り返してみましょう。
Posted in UT

広域システム 09.28

x=rnorm(n=1, mean=100, sd=15)
y=runif(n=3, min=1, max=10)
N = 10000

random.data = rnorm(N, mean=0, sd=1)
hist(random.data, nclass = 50, col = "navy", xlab = "Data",
     probability = T, main = "Histogram of Random Data")
dens = density(random.data)
lines(dens, col = "orange", lwd = 4)

sample(1:10,3)
sample(c(“gu”,“choki”,“pa”),1)
sample(1:10)
sample(0:1, 10, replace=T)

# FOR loop
for (i_loop in 1:5) {
  print(i_loop)
}

for (i_loop in 1:5) {
  print(c(i_loop, 2^i_loop))
}

counter <- 1
while(counter<=10){
  print(counter)
  counter<-counter+1
}

counter <- 1
while(counter^2 <= 10){
  print(c(counter, counter^2))
  counter<-counter+1
}

affil<-"cogsci"
if (affil=="cogsci") {
  print("you are wonderful")
}

affil<-"phil"
if (affil=="cogsci") {
  print("you are wonderful")
} else {
  print("still, you are wonderful")
}

counter=6
repeat{
  print(counter)
  counter = counter + 1
  if(counter>5){break}
}

counter=6
repeat{
  if(counter>5){break}
  print(counter)
  counter+counter+1
}

six.counter=0; N = 1000
for (i_loop in 1:N) {
  die<-sample(1:6,1)
  if (die==6) {six.counter=six.counter+1}
}
six.counter/N

N=1000; six.counter=rep(0,N);
for (i_loop in 1:N) {
  die<-sample(1:6,1)
  if (die==6) {six.counter[i_loop]=1}
}
plot(1:N,cumsum(six.counter)/(1:1000),type='l',ylim=c(0,1),
     lwd=2)
abline(h=1/6,lwd=2,col='red')

N = 1000
die.all <- sample(1:6,N,replace=T)
six.index <- die.all==6
par(mfrow = c(2,1))
par(oma=c(2,2,0,0),mar=c(4,4,1,1),mfrow=c(2,1))
plot(1:N, die.all, pch=20, col = 'red', ylim = c(0,7),
     ylab = "Result", xlab = "trial")
plot(1:N,cumsum(six.index)/(1:1000), type='l', ylim=c(0,1), lwd=2,
     ylab = "P(die = 6)", xlab = "trial")
abline(h=1/6,lwd=2,col='red')

# CLT w/ loop
N=10;nRep=10000;
means<-rep(0,nRep)
for (i_rep in 1:nRep) {
  dat<-runif(N)
  means[i_rep]=mean(dat)
}
hist(means,nclass=50,probability=T)
dens<-density(means)
lines(dens,col='skyblue',lwd=3)
xs=seq(-0,1,0.01)
theo.dens<-dnorm(xs,mean=0.5,sd=sqrt((1/12)/N))
lines(xs,theo.dens,col='orange',lwd=3,lty=2)

# CLT w/o loop
N=10
nRep=10000
dat<-matrix(runif(N*nRep),nrow=N)
means<-colMeans(dat)
hist(means,nclass=50,probability=T)
dens<-density(means); lines(dens,col='skyblue',lwd=3)
xs=seq(-0,1,0.01)
theo.dens<-dnorm(xs,mean=0.5,sd=sqrt((1/12)/N))
lines(xs,theo.dens,col='orange',lwd=3,lty=2)

r=-99;v1=1633;v2=355
while (r!=0){
  r=v1%%v2
  print(paste('v1 =',v1,', v2 = ',v2,',remainder = ',r))
  v1=v2
  v2=r
}

tol = 1e-7;grad = 1e10; lambda = 0.1;x = 10; x.hist = x
grad = 2*x+2
while (abs(grad)>tol){
  x = x - lambda*grad
  x.hist=c(x.hist,x)
  grad = 2*x+2
}
x.temp=seq(-10,10,0.1)
plot(x.temp, x.temp^2+2*x.temp,type='l',lwd=2)
lines(x.hist,x.hist^2+2*x.hist,type='o',pch=1,col='red',cex=1)
Posted in UT

データ解析基礎論a 分散分析3

source("http://www.matsuka.info/univ/course_folder/cutil.R")

dat<-read.csv("http://www.matsuka.info/data_folder/dktb316.txt")
dat.aov<-aov(words~method+subj+Error(subj:method),data=dat)
dat.aov2<-aov(words~method+Error(subj+subj:method),data=dat)
dat.aov.BTW <-aov(words~method,data=dat)

summary(dat.aov)
RB.qtukey(dat.aov,dat, 0.05)

dat<-read.csv("http://www.matsuka.info/data_folder/dktb3211.txt")
interaction.plot(dat$duration,  # x軸
                 dat$method,    # まとめる変数    
                 dat$result,    # y軸 
                 pch=c(20,20), 
                 col=c("skyblue","orange"),
                 ylab="score",
                 xlab="Duration",
                 lwd=3, 
                 type='b',
                 cex=2, 
                 trace.label="Method")

dat.aov <- aov(result~method*duration+Error(s+s:duration),dat)
TSME<-SPF.tsme(dat.aov,dat,"result")

dat<-read.csv("http://www.matsuka.info/data_folder/dktb3218.txt")
dat.aov<-aov(result~method+duration+method:duration + Error(s+method:s+duration:s+method:duration:s), data=dat)
TSME<-RBF.tsme(dat.aov, dat, "result")

認知情報解析演習 居住区問題

n.circle = 20; n.sharp = 20; size = 10
loc = sample(1:size^2, n.circle+n.sharp)
type = c(rep(1,n.circle),rep(2,n.sharp))
# circle = 1; sharp = 2
people = cbind(type,loc)
state.mat = matrix(0,size,size)
state.mat[people[,2]]=people[,1]

p.move = #1/(2*n.circle)
p.other = 0.1
c.count = 3

idx = cbind(rep(1:size,size),sort(rep(1:size,size)))

#
term = F
while (term == F){
  term = T
  rand.order = sample(1:nrow(people))
  for (i.p in 1:nrow(people)){
    if (people[rand.order[i.p],1]==1){
      # circle
      if (runif(1) < p.move){
        empty = 1:(size^2)
        empty = empty[-sort(people[,2])]
        people[rand.order[i.p],2] = sample(empty,1)
        state.mat = matrix(0,size,size)
        state.mat[people[,2]]=people[,1]
        term = F
      } 
    } else {
      # sharp
      x.min = max(idx[people[rand.order[i.p],2],1]-1,1)
      x.max = min(idx[people[rand.order[i.p],2],1]+1,size)
      y.min = max(idx[people[rand.order[i.p],2],2]-1,1)
      y.max = min(idx[people[rand.order[i.p],2],2]+1,size)
      circle.in = sum(state.mat[x.min:x.max,y.min:y.max]==1)
      if (circle.in >= c.count){
        empty = 1:(size^2)
        empty = empty[-sort(people[,2])]
        people[rand.order[i.p],2] = sample(empty,1)
        state.mat = matrix(0,size,size)
        state.mat[people[,2]]=people[,1]
        term = F
        #print('moved')
      }
    }
  }
  for (i.p in 1:nrow(people)){
    
    if (people[rand.order[i.p],1]==2){
      x.min = max(idx[people[i.p,2],1]-1,1)
      x.max = min(idx[people[i.p,2],1]+1,size)
      y.min = max(idx[people[i.p,2],2]-1,1)
      y.max = min(idx[people[i.p,2],2]+1,size)
      circle.in = sum(state.mat[x.min:x.max,y.min:y.max]==1)
      print(circle.in)
      if (circle.in >= c.count){
        term = F
        break
      }
    }
  }
}
plot(0,0, type= 'n', xlim = c(0,11),ylim=c(0,11))
lab = c("0","#")
text(idx[people[,2],1],idx[people[,2],2],lab[people[,1]],cex=3)
ab = seq(0.5,10.5,1)
for (i in 1:11){
  abline(h=ab[i],col='red')
  abline(v=ab[i],col='red')
}

認知情報解析演習 石とりゲーム(bugあり)

col1 = matrix(c(rep(0,4),c(1,0,0,0),c(1,1,0,0),c(1,1,1,0),rep(1,4)),nrow=4,byrow=T)
col2 = matrix(c(rep(10,4),c(11,10,10,10),c(11,11,10,10),c(11,11,11,10),rep(11,4)),nrow=4,byrow=T)
col3 = matrix(c(rep(100,4),c(101,100,100,100),c(101,101,100,100),c(101,101,101,100),rep(101,4)),nrow=4,byrow=T)
act.list = list()
state.temp = list()
counter = 0
Q1 = list()
Q2 = list()
for (i.c1 in 1:5){
  if (sum(col1[,i.c1])==0){
    act1 = c()
  } else {
    act1 = seq(1,sum(col1[,i.c1]),1)
  }
  for (i.c2 in 1:5){
    if (sum(col2[,i.c2])==40){
      act2 = c()
    } else {
      act2 = seq(11,sum(col2[,i.c2]==11)*11,11)
    }
    for (i.c3 in 1:5){
      if (sum(col3[,i.c3])==400){
        act3 = c()
      } else {
        act3 = seq(101,sum(col3[,i.c3]==101)*101,101)
      }
      counter = counter + 1
      state.temp[[counter]] = cbind(col1[,i.c1],col2[,i.c2],col3[,i.c3])
      act.list[[counter]] = c(act1,act2,act3)
      Q1[[counter]] = rep(0, length(c(act1,act2,act3)))
      Q2[[counter]] = rep(0, length(c(act1,act2,act3)))
    }
  }
}
rm.stone <- function(act, st.shape){
  if (act == -99){s}
  if (act > 100){
    n.remove = act%%100
    n.stone = length(which(st.shape[,3]==101))
    start = (5-n.stone)
    st.shape[(start:(start+n.remove-1)),3] = 100
  } else {
    if (act > 10){
      n.remove = act%%10
      n.stone = length(which(st.shape[,2]==11))
      start = (5-n.stone)
      st.shape[(start:(start+n.remove-1)),2] = 10
    } else {
      n.remove = act
      n.stone = length(which(st.shape[,1]==1))
      start = (5-n.stone)
      st.shape[(start:(start+n.remove-1)),1] = 0
    }
  }
  return(st.shape)
}

id.state <- function(st.shape, state.temp){
  for (i.st in 1:125){
    if  (all(st.shape == state.temp[[i.st]])){
      state.idx = i.st 
      break
    }
  }
  return(state.idx)
}

ck.act <- function(Q, act.vec, eta){
  if (is.null(act.vec)){
    return(list(act = -99, act.idx = -99))
    break
  }
  if (length(act.vec)==1){
    act = act.vec
  } else {
    p = exp(Q[[state.idx]])/sum(exp(Q[[state.idx]]))
    act = sample(act.vec, 1, prob = p)
  }
  act.idx = which(act.vec==act)
  return(list(act = act, act.idx = act.idx))
}


gamma=1;alpha = 0.1;n.rep=10000
for (i.rep in 1:n.rep){
  # first action
  state.idx = 125; counter = 1
  st.shape = state.temp[[state.idx]]
  res.act = ck.act(Q1,act.list[[state.idx]],eta)
  act = res.act$act;act.idx = res.act$act.idx
  state.old = state.idx
  act.old = act.idx
  
  # 2nd to last
  while (state.idx != 1) {
    counter = counter + 1
    st.shape <- rm.stone(act, st.shape)
    state.idx <- id.state(st.shape, state.temp)
    if (counter%%2==1) {
      res.act = ck.act(Q1,act.list[[state.idx]],eta)
    } else {
      res.act = ck.act(Q2,act.list[[state.idx]],eta)
    }
    act = res.act$act; act.idx = res.act$act.idx
    if (state.idx == 1){
      if (counter%%2==1){rew1 = -1; rew2 = 1;} else {rew1 = 1; rew2 = -1;}
      Q1[[state.old]][act.old]=Q1[[state.old]][act.old]
         +alpha*(rew1-Q1[[state.old]][act.old])
      Q2[[state.old]][act.old]=Q2[[state.old]][act.old]
         +alpha*rew2-Q2[[state.old]][act.old])
    } else {
      rew1 = 0; 
      rew2 =0;
      if (counter%%2==1){
        Q1[[state.old]][act.old]=Q1[[state.old]][act.old]
          +alpha*(rew1+gamma* Q1[[state.idx]][act.idx]-Q1[[state.old]][act.old])
      } else {
        Q2[[state.old]][act.old]=Q2[[state.old]][act.old]
          +alpha*(rew2+gamma* Q2[[state.idx]][act.idx]-Q2[[state.old]][act.old])
      }
    }
    
    state.old = state.idx
    act.old = act.idx
  }
}  

2019年度 データ解析基礎論a 分散分析2

source("http://www.matsuka.info/univ/course_folder/cutil.R")
adj.alpha(5,0.05)

# anova
dat<-read.csv("http://www.matsuka.info/data_folder/dktb312.csv")
dat$method=factor(dat$method, levels(dat$method)[c(1,3,4,2)])
dat.aov<-aov(result~method, data=dat)
summary(dat.aov)

# multiple comparison 
# do not use these command
# use simpler one  - "cu.bonF1F"
dat.means<-tapply(dat$result,dat$method,mean)
new.alpha = 1-(1-0.05)^(1/5)
cont=c(-3,1,1,1)
bunshi=sum(cont*dat.means)
bunbo=sqrt(5.29*(sum((cont^2)/8)))
t.value=bunshi/bunbo
2*(1-pt(t.value,28))


# cu.bonF1F
new.alpha<-adj.alpha(5,0.05) 
cu.bonF1F(dat.aov,dat,c(-3,1,1,1),new.alpha)
cu.bonF1F(dat.aov,dat,c(-1,1,0,0),new.alpha)
cu.bonF1F(dat.aov,dat,c(-1,0,1,0),new.alpha)
cu.bonF1F(dat.aov,dat,c(-1,0,0,1),new.alpha)
cu.bonF1F(dat.aov,dat,c(0,-2,1,1),new.alpha)

# cu.scheffe1F
new.f<-adj.f(4,3,28,0.05)
cu.scheffe1F(dat.aov,dat,c(-3,1,1,1))


# 2 way ANOVA 
dat <- read.csv("http://peach.l.chiba-u.ac.jp/course_folder/datW03R.txt")
interaction.plot(dat$gender,
                 dat$affil,
                 dat$shoesize, 
                 pch=c(20,20), 
                 col=c("skyblue","orange"), 
                 xlab="gender", ylab="shoesize", 
                 lwd=3,type='b',cex=2,
                 trace.label="Affiliation")
dat.aov=aov(shoesize~gender*affil, data=dat)
dat.aov.sum=summary(dat.aov)

# testing simple main effect
# do not use these command
# use simpler one - CRF.tsme
means<-tapply(dat$shoesize, list(dat$gender,dat$affil), mean)
SS_gen_CS<- 5*(means[2,1]^2 + means[1,1]^2 -0.5*sum(means[,1])^2) # SS_gender CS
SS_gen_PS<- 5*(means[2,2]^2 + means[1,2]^2 -0.5*sum(means[,2])^2) # SS_gender PS
dat.aov.sum=summary(dat.aov)   # ANOVA table
MSe=dat.aov.sum[[1]][4,3]      # MSE from ANOVA table or MSe=0.62
dfE=dat.aov.sum[[1]][4,1]      # DF for error or dfE=16
dfG=1                          # DF for gender
F_gen_CS=(SS_gen_CS/dfG)/MSe   # F-value for gender effect given CS
F_gen_PS=(SS_gen_PS/dfG)/MSe   # F-value for gender effect given PS
P_gen_CS=1-pf(F_gen_CS,1,dfE)  # p-value for gender effect given CS
P_gen_PS=1-pf(F_gen_PS,1,dfE)  
SS_affil_F<- 5*(means[1,1]^2+means[1,2]^2-0.5*sum(means[1,])^2) #SS_affil | F
SS_affil_M<- 5*(means[2,1]^2+means[2,2]^2-0.5*sum(means[2,])^2) #SS_affil | M
dfA=1				          # DF for affil
F_affil_F=SS_affil_F/dfA/MSe         # F-value for affiliation effect | F
F_affil_M=SS_affil_M/dfA/MSe         # F-value for affiliation effect | M
P_affil_F=1-pf(F_affil_F,1,dfE)      # p-value for affiliation effect | F
P_affil_M=1-pf(F_affil_M,1,dfE)      # p-value for affiliation effect | M

# testint simple main effect w/ CRF.tsme
tsme = CRF.tsme(dat.aov, dat)
             
# another 2-way ANOVA
dat <-read.csv("http://www.matsuka.info/data_folder/dktb321.txt")
interaction.plot(dat$duration,dat$method,dat$result, pch=c(20,20), 
                 col=c("skyblue","orange"), ylab="score", xlab="Duration",
                 lwd=3,type='b',cex=2,trace.label="Method")
mod1=aov(result~method+duration,data=dat)
mod1.sum=print(summary(mod1))

mod2=aov(result~method*duration,data=dat)
mod2.sum=print(summary(mod2))
CRF.tsme(mod2, dat)


# plotting 
dat<-read.csv("http://www.matsuka.info/univ/course_folder/datW03R.txt",header=T)
means<-tapply(dat$shoesize,list(dat$gender, dat$affil),mean)
Ns<-tapply(dat$shoesize,list(dat$gender, dat$affil),length)
sds<-tapply(dat$shoesize,list(dat$gender, dat$affil),sd)
sems<-sds/sqrt(Ns)
plot(c(0,1),means[,1],type='o',col='skyblue', xlim=c(-0.1,1.1), lwd=2, cex=2, pch=20,
     ylim=c(min(means)*0.975, max(means)*1.025), xlab="gender", ylab="shoesize", xaxt="n")
axis(1,c(0,1),c("female","male"))
lines(c(0,1),means[,2],type='o',col='orange',lwd=2,cex=2,pch=20)
lines(c(0,0),c(means[1,1]-sems[1,1],means[1,1]+sems[1,1]),col="skyblue",lwd=2)
lines(c(0,0),c(means[1,2]-sems[1,2],means[1,2]+sems[1,2]),col="orange",lwd=2)
lines(c(1,1),c(means[2,1]-sems[2,1],means[2,1]+sems[2,1]),col="skyblue",lwd=2)
lines(c(1,1),c(means[2,2]-sems[2,2],means[2,2]+sems[2,2]),col="orange",lwd=2)
legend("topleft",c("CogSci","PsySci"),col=c("skyblue","orange"),lwd=2)

院:認識情報解析

library(rjags)
source("http://peach.l.chiba-u.ac.jp/course_folder/HDI_revised.txt")

dat<-read.csv("http://www.matsuka.info/data_folder/HtWtData110.csv")
library(plot3D)
w = seq(80,360,length.out=100)
h = seq(50, 75, length.out=100)
M <- mesh(w,h)
P.male = 1/(1+exp(-1*(0.018*M$x+0.7*M$y-50)))

scatter3D(dat$weight, dat$height, dat$mal, pch = 19, cex = 2,
          theta = 30, phi = 45, ticktype = "detailed", zlim=c(-0.1,1),ylim=c(50,78),xlim=c(80,360),
          xlab = "weight", ylab = "height", zlab = "P(male)",
          surf = list(x = M$x, y = M$y, z = P.male,facets = NA))


y = dat$male; x = dat$weight; Ntotal = length(y)
dataList = list(y = y, x = x, Ntotal = Ntotal)

model.txt = "
data {
  xm <- mean(x)
  xsd <- sd(x)
  for (i in 1:Ntotal){
    zx[i] = (x[i] - xm)/xsd
  }
}
model {
  for ( i_data in 1:Ntotal ) {
    y[ i_data ] ~ dbern(ilogit( zbeta0 + zbeta * zx[i_data]))
  }
  zbeta0 ~ dnorm(0, 1/2^2)
  zbeta ~ dnorm(0, 1/2^2)

  beta <- zbeta / xsd
  beta0 <- zbeta0 - zbeta * xm/xsd
}"
writeLines(model.txt, "model1.txt")
parameters = c( "beta0" ,  "beta")
jagsModel = jags.model( "model1.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$weight,dat$male,xlim=c(90,280),yaxt="n",ylab="Male / Female",
     xlab="Weight", cex=2.5)
axis(2,at = 0:1,labels=c("Femal","Male"))
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
temp.x = seq(90,280,length.out = 100)
for (i_sample in 1:n2plot) {
  temp.y = 1/(1+exp(-1*(mcmcMat[idx[i_sample],2] + mcmcMat[idx[i_sample],1]*temp.x)))
  lines(temp.x, temp.y, col='orange', lwd=2)
}


x = cbind(dat$weight,dat$height);Nx = ncol(x)
dataList = list(y = y, x = x, Ntotal = Ntotal, Nx = Nx)

model.txt = "
data {
  for (j in 1:Nx){
    xm[j] <- mean(x[,j])
    xsd[j] <- sd(x[,j])
    for (i in 1:Ntotal){
      zx[i,j] = (x[i,j] - xm[j])/xsd[j]
    }
  }
}
model {
  for ( i_data in 1:Ntotal ) {
    y[ i_data ] ~ dbern(ilogit( zbeta0 + sum(zbeta[1:Nx] * zx[i_data, 1:Nx ])))
  }
  zbeta0 ~ dnorm(0, 1/2^2)
  for (j in 1:Nx){
    zbeta[j] ~ dnorm(0, 1/2^2)
  }
  beta[1:Nx] <- zbeta[1:Nx] / xsd[1:Nx]
  beta0 <- zbeta0 -sum(zbeta[1:Nx] * xm[1:Nx]/xsd[1:Nx])
}"
writeLines(model.txt, "model.txt")

parameters = c( "beta0" ,  "beta")
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)
par(mfrow=c(1,3))
HDI.plot(mcmcMat[,3],xlabel='intercept')
HDI.plot(mcmcMat[,1],xlabel='weight')
HDI.plot(mcmcMat[,2],xlabel='height')

par(mfrow=c(1,1))
plot(dat$weight,dat$height,xlab="Weight", ylab="Height", type="n")
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
for (i_sample in 1:n2plot) {
  abline(a=-1*mcmcMat[idx[i_sample],3]/mcmcMat[idx[i_sample],2],
         b=-1*mcmcMat[idx[i_sample],1]/mcmcMat[idx[i_sample],2],col="orange")

}
points(dat$weight,dat$height,pch=paste(dat$male), cex=1.5)

# un-even data
x = rnorm(300)
pr = 1/(1+exp(2*x))
y = pr < runif(300)
plot(x,y)

remove.id = sample(which(y == 0),120)

Ntotal = length(y[-remove.id])
dataList = list(y = y[-remove.id], x = x[-remove.id], Ntotal = Ntotal)

model.txt = "
data {
xm <- mean(x)
xsd <- sd(x)
for (i in 1:Ntotal){
zx[i] = (x[i] - xm)/xsd
}
}
model {
for ( i_data in 1:Ntotal ) {
y[ i_data ] ~ dbern(ilogit( zbeta0 + zbeta * zx[i_data]))
}
zbeta0 ~ dnorm(0, 1/2^2)
zbeta ~ dnorm(0, 1/2^2)

beta <- zbeta / xsd
beta0 <- zbeta0 - zbeta * xm/xsd
}"
writeLines(model.txt, "model1.txt")
parameters = c( "beta0" ,  "beta")
jagsModel = jags.model( "model1.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(x[-remove.id],y[-remove.id],xlim=c(-3,3))
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
temp.x = seq(-3,3,length.out = 100)
for (i_sample in 1:n2plot) {
  temp.y = 1/(1+exp(-1*(mcmcMat[idx[i_sample],2] + mcmcMat[idx[i_sample],1]*temp.x)))
  lines(temp.x, temp.y, col='orange', lwd=2)
}


x1 = rnorm(150)
x2 = x1*0.9+rnorm(150,0,0.5)
pr = 1/(1+exp(x1+x2))
y = pr < runif(150)
Ntotal = length(y)
dataList = list(y = y, x = cbind(x1,x2), Ntotal = Ntotal, Nx = 2)


model.txt = "
data {
for (j in 1:Nx){
xm[j] <- mean(x[,j])
xsd[j] <- sd(x[,j])
for (i in 1:Ntotal){
zx[i,j] = (x[i,j] - xm[j])/xsd[j]
}
}
}
model {
for ( i_data in 1:Ntotal ) {
y[ i_data ] ~ dbern(ilogit( zbeta0 + sum(zbeta[1:Nx] * zx[i_data, 1:Nx ])))
}
zbeta0 ~ dnorm(0, 1/2^2)
for (j in 1:Nx){
zbeta[j] ~ dnorm(0, 1/2^2)
}
beta[1:Nx] <- zbeta[1:Nx] / xsd[1:Nx]
beta0 <- zbeta0 -sum(zbeta[1:Nx] * xm[1:Nx]/xsd[1:Nx])
}"
writeLines(model.txt, "model.txt")

parameters = c( "beta0" ,  "beta")
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)
plot(x1,x2,xlab="x1", ylab="x2", type="n")
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
for (i_sample in 1:n2plot) {
  abline(a=-1*mcmcMat[idx[i_sample],3]/mcmcMat[idx[i_sample],2],
         b=-1*mcmcMat[idx[i_sample],1]/mcmcMat[idx[i_sample],2],col="orange")

}
points(x1,x2,pch=paste(y), cex=1.5)

# guessing
y = dat$male; x = dat$weight; Ntotal = length(y)
dataList = list(y = y, x = x, Ntotal = Ntotal)

model.txt = "
data {
xm <- mean(x)
xsd <- sd(x)
for (i in 1:Ntotal){
zx[i] = (x[i] - xm)/xsd
}
}
model {
for ( i_data in 1:Ntotal ) {
y[ i_data ] ~ dbern(mu[i_data])
mu[i_data] <- (guess*0.5 + (1-guess)*ilogit( zbeta0 + zbeta * zx[i_data]))
}
zbeta0 ~ dnorm(0, 1/2^2)
zbeta ~ dnorm(0, 1/2^2)
guess ~ dbeta(1,9)
beta <- zbeta / xsd
beta0 <- zbeta0 - zbeta * xm/xsd
}"
writeLines(model.txt, "model1.txt")
parameters = c( "beta0" ,  "beta", "guess")
jagsModel = jags.model( "model1.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(x[-remove.id],y[-remove.id],xlim=c(-3,3))
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
temp.x = seq(-3,3,length.out = 100)
for (i_sample in 1:n2plot) {
  temp.y = 1/(1+exp(-1*(mcmcMat[idx[i_sample],2] + mcmcMat[idx[i_sample],1]*temp.x)))
  lines(temp.x, temp.y, col='orange', lwd=2)
}
par(mfrow=c(1,3))
HDI.plot(mcmcMat[,2],xlabel='intercept')
HDI.plot(mcmcMat[,1],xlabel='weight')
HDI.plot(mcmcMat[,3],xlabel='guessing')

par(mfrow=c(1,1))
plot(dat$weight,dat$male,xlim=c(90,280),yaxt="n",ylab="Male / Female",
     xlab="Weight", cex=2.5)
axis(2,at = 0:1,labels=c("Femal","Male"))
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
temp.x = seq(90,280,length.out = 100)
for (i_sample in 1:n2plot) {
  temp.y = mcmcMat[idx[i_sample],3]/2+(1-mcmcMat[idx[i_sample],3])*1/(1+exp(-1*(mcmcMat[idx[i_sample],2] + mcmcMat[idx[i_sample],1]*temp.x)))
  lines(temp.x, temp.y, col='orange', lwd=2)
}

# nomial predictors
model.txt = "
model {
  for ( i.data in 1:Ntotal ) {
    y[ i.data ] ~ dbin(mu[i.data],N[i.data])
    mu[i.data] ~ dbeta(omega[x[i.data]]*(kappa-2)+1,(1-omega[x[i.data]])*(kappa-2)+1)
  }
  for (i.pos in 1:Npos){
    omega[i.pos] <- ilogit(a0+a[i.pos])
    a[i.pos] ~ dnorm(0.0, 1/aSigma^2)
  }
  a0 ~  dnorm(0,1/2^2)
  aSigma ~ dgamma(1.64, 0.32)
  kappa <- kappaMinusTwo +2
  kappaMinusTwo ~ dgamma(0.01,0.01)
  for (i.pos in 1:Npos){
    m[i.pos] <- a0+a[i.pos]
  }
  b0 <- mean(m[1:Npos])
  for (i.pos in 1:Npos){
    b[i.pos] <- m[i.pos] - b0
  }
}"
writeLines(model.txt, "model.txt")

dat<-read.csv("http://www.matsuka.info/data_folder/BattingAverage.csv")
y = dat$Hits
N = dat$AtBats
x = dat$PriPosNumber
Ntotal = length(y)
Npos = length(unique(x))
dataList = list(y = y, x = x, N = N, Ntotal = Ntotal, Npos = Npos)
parameters = c( "b0" ,  "b", "omega")
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)

par(mfrow=c(3,3))
for (i.pos in 1:9){
  HDI.plot(mcmcMat[,i.pos+10])
}

par(mfrow=c(2,2))
HDI.plot(mcmcMat[,1]-mcmcMat[,2])
HDI.plot(mcmcMat[,2]-mcmcMat[,3])
HDI.plot(mcmcMat[,11]-mcmcMat[,12])
HDI.plot(mcmcMat[,12]-mcmcMat[,13])

# softmax regression
x1 = runif(500, min=-2, max = 2)
x2 = runif(500, min=-2, max = 2)
b0 = c(0,-3,-4,-5)
b1 = c(0,-5,-1,10)
b2 = c(0,-5,10,-1)
l1 = b0[1]+b1[1]*x1+b2[1]*x2
l2 = b0[2]+b1[2]*x1+b2[2]*x2
l3 = b0[3]+b1[3]*x1+b2[3]*x2
l4 = b0[4]+b1[4]*x1+b2[4]*x2
p1 = exp(l1)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p2 = exp(l2)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p3 = exp(l3)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p4 = exp(l4)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
ps = cbind(p1,p2,p3,p4)
y = apply(ps,1,which.max)
plot(x1,x2,pch=y,col=y)

b0 = c(0,-4,-1,-1)
b1 = c(0,-5,1,3)
b2 = c(0,0,-5,3)
l1 = b0[1]+b1[1]*x1+b2[1]*x2
l2 = b0[2]+b1[2]*x1+b2[2]*x2
l3 = b0[3]+b1[3]*x1+b2[3]*x2
l4 = b0[4]+b1[4]*x1+b2[4]*x2
p1 = exp(l1)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p2 = exp(l2)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p3 = exp(l3)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p4 = exp(l4)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
ps = cbind(p1,p2,p3,p4)
y = apply(ps,1,which.max)
plot(x1,x2,pch=y,col=y)

p1 = exp(l1)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p2 = exp(l2)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p3 = exp(l3)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p4 = exp(l4)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
ps = cbind(p1,p2,p3,p4)
p12 = pmax(p1,p2)
p34 = pmax(p3,p4)
y12vs34 = apply(cbind(p1,p2),1,which.max)
plot(x1,x2,pch=y12vs34,col=y12vs34)
y1vs2 = apply(cbind(p1,p3),1,which.max)
points(x1,x2,pch=y1vs2+2,col=y1vs2+2)
y3vs4 = apply(cbind(p1,p4),1,which.max)
points(x1,x2,pch=y3vs4+6,col=y3vs4+6)



model.txt = "
data {
  for ( j in 1:Nx ) {
    xm[j]  <- mean(x[,j])
    xsd[j] <-   sd(x[,j])
    for ( i in 1:Ntotal ) {
      zx[i,j] <- ( x[i,j] - xm[j] ) / xsd[j]
    }
  }
}
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dcat(mu[1:Nout,i])
    mu[1:Nout,i] <- explambda[1:Nout,i]/sum(explambda[1:Nout,i])
    for (k in 1:Nout){
      explambda[k,i]=exp(zbeta0[k] + sum(zbeta[k,1:Nx] * zx[i, 1:Nx ]))
    }
  }
  zbeta0[1] = 0
  for (j in 1:Nx){
    zbeta[1,j] <- 0
  }
  for (k in 2:Nout){
    zbeta0[k] ~ dnorm(0, 1/2^2)
    for (j in 1:Nx){
      zbeta[k,j]~dnorm(0, 1/2^2)
    }
  }
  for ( k in 1:Nout ) {
    beta[k,1:Nx] <- zbeta[k,1:Nx] / xsd[1:Nx]
    beta0[k] <- zbeta0[k] - sum( zbeta[k,1:Nx] * xm[1:Nx] / xsd[1:Nx] )
  }
}"
writeLines(model.txt, "model.txt")

dat<-read.csv( "http://www.matsuka.info/data_folder/CondLogistRegData1.csv" )
y = dat$Y
x = cbind(dat[,1],dat[,2])
Ntotal = length(y)
Nout = length(unique(y))
dataList = list(y = y, x = x, Nx = 2, Ntotal = Ntotal, Nout = Nout)
parameters = c( "beta0" ,  "beta")
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)
par(mfrow=c(1,3))
HDI.plot(mcmcMat[,7+0],xlab='intercept')
HDI.plot(mcmcMat[,1+0],xlab='b1')
HDI.plot(mcmcMat[,4+0],xlab='b2')

model = "
data {
  for ( j in 1:Nx ) {
    xm[j]  <- mean(x[,j])
    xsd[j] <-   sd(x[,j])
    for ( i in 1:Ntotal ) {
      zx[i,j] <- ( x[i,j] - xm[j] ) / xsd[j]
    }
  }
}
# Specify the model for standardized data:
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dcat( mu[1:Nout,i] )
    mu[1,i] <- phi[1,i]
    mu[2,i] <- phi[2,i] * (1-phi[1,i])
    mu[3,i] <- phi[3,i] * (1-phi[2,i]) * (1-phi[1,i])
    mu[4,i] <- (1-phi[3,i]) * (1-phi[2,i]) * (1-phi[1,i])
    for ( r in 1:(Nout-1) ) {
      phi[r,i] <- ilogit( zbeta0[r] + sum( zbeta[r,1:Nx] * zx[i,1:Nx] ) )
    }
  }
  for ( r in 1:(Nout-1) ) {
    zbeta0[r] ~ dnorm( 0 , 1/20^2 )
    for ( j in 1:Nx ) {
      zbeta[r,j] ~ dnorm( 0 , 1/20^2 )
    }
  }
  for ( r in 1:(Nout-1) ) {
    beta[r,1:Nx] <- zbeta[r,1:Nx] / xsd[1:Nx]
    beta0[r] <- zbeta0[r] - sum( zbeta[r,1:Nx] * xm[1:Nx] / xsd[1:Nx] )
  }
}
"
writeLines( model , "model.txt" )

dat<-read.csv( "http://www.matsuka.info/data_folder/SoftmaxRegData2.csv" )
y = dat$Y
x = cbind(dat[,1],dat[,2])
Ntotal = length(y)
Nout = length(unique(y))
dataList = list(y = y, x = x, Nx = 2, Ntotal = Ntotal, Nout = Nout)
parameters = c( "beta0" ,  "beta")
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)
par(mfrow=c(1,1))
plot(x[,1],x[,2],col=y)
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
temp.x = seq(-3,3,length.out = 100)
for (i.cat in 0:2){
  for (i_sample in 1:n2plot) {
    abline(a=-1*mcmcMat[idx[i_sample],7+i.cat]/mcmcMat[idx[i_sample],4+i.cat],
           b=-1*mcmcMat[idx[i_sample],1+i.cat]/mcmcMat[idx[i_sample],4+i.cat],col="orange")
  }
}
model2 = "
data {
  for ( j in 1:Nx ) {
    xm[j]  <- mean(x[,j])
    xsd[j] <-   sd(x[,j])
      for ( i in 1:Ntotal ) {
        zx[i,j] <- ( x[i,j] - xm[j] ) / xsd[j]
      }
    }
  }
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dcat( mu[1:Nout,i] )
    mu[1,i] <- phi[2,i] * phi[1,i]
    mu[2,i] <- (1-phi[2,i]) * phi[1,i]
    mu[3,i] <- phi[3,i] * (1-phi[1,i])
    mu[4,i] <- (1-phi[3,i]) * (1-phi[1,i])
    for ( r in 1:(Nout-1) ) {
      phi[r,i] <- ilogit( zbeta0[r] + sum( zbeta[r,1:Nx] * zx[i,1:Nx] ) )
    }
  }
  for ( r in 1:(Nout-1) ) {
    zbeta0[r] ~ dnorm( 0 , 1/20^2 )
    for ( j in 1:Nx ) {
      zbeta[r,j] ~ dnorm( 0 , 1/20^2 )
    }
  }
  for ( r in 1:(Nout-1) ) {
    beta[r,1:Nx] <- zbeta[r,1:Nx] / xsd[1:Nx]
    beta0[r] <- zbeta0[r] - sum( zbeta[r,1:Nx] * xm[1:Nx] / xsd[1:Nx] )
  }
}"
writeLines( modelString , con="TEMPmodel.txt" )