# 認知情報解析学 MCMC

```import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import seaborn as sb

# one variable
N = 20; z = 14
a = 1; b = 1
sigma = 0.2
nIter = 100000
theta = np.random.rand(1)
thetaHist = np.zeros(nIter)
for i_iter in range(nIter):
proposed = theta + np.random.normal(0,sigma,1)
while proposed > 1 or proposed < 0:
proposed = theta + np.random.normal(1)*sigma
vProp = (proposed**z * (1-proposed)**(N-z) * proposed**(a-1) * (1-proposed)**(b-1) )
vCurr = (theta**z * (1-theta)**(N-z) * theta**(a-1) * (1-theta)**(b-1) )
pmove = min(1,vProp/vCurr)
if np.random.rand(1) < pmove:
theta = proposed
thetaHist[i_iter] = theta
sb.distplot( thetaHist, bins=100 )

# two variable
N1 = 8; z1 = 6
N2 = 7; z2 = 2
a1 = 2; b1 = 2
a2 = 2; b2 = 2

sigma = 0.2
nIter = 1000000
theta1 = np.random.rand(1)
theta2 = np.random.rand(1)
thetaHist = np.zeros([nIter,2])
for i_iter in range(nIter):

prop1 = theta1 + np.random.normal(0,sigma,1)
while prop1 > 1 or prop1 < 0:
prop1 = theta1 + np.random.normal(0,sigma,1)

prop2 = theta2 + np.random.normal(0,sigma,1)
while prop2 > 1 or prop2 < 0:
prop2 = theta2 + np.random.normal(0,sigma,1)

vProp = (prop1**z1* (1-prop1)**(N1-z1)* prop2**z2 *(1-prop2)**(N2-z2)
*prop1**(a1-1) *(1-prop1)**(b1-1) *prop2**(a2-1) *(1-prop2)**(b2-1))
vCurr = (theta1**z1*(1-theta1)**(N1-z1)*theta2**z2*(1-theta2)**(N2-z2)
*theta1**(a1-1)*(1-theta1)**(b1-1)*theta2**(a2-1)*(1-theta2)**(b2-1))
pmove = min(1,vProp/vCurr)
if np.random.rand(1) < pmove:
theta1 = prop1
theta2 = prop2
thetaHist[i_iter,...] = [theta1,theta2]
print(np.mean(thetaHist[...,1]),np.mean(thetaHist[...,0]))
sb.jointplot(x=thetaHist[...,0], y=thetaHist[...,1], kind='kde',space=0, ratio=4)
```

# iterative prisoners’ dilemma

```pmat = np.array([[-25,0],[-20,-5]])

def iterPD(action1, action2, pmat):
# A = defect & defect (-10); B = defect & coop (1)
# C = coop & defect   (9)  ;   D = coop & coop (11)

result = action1 + action2
if result == -10:
po = [pmat[0,0],pmat[0,0]]
elif result == -9:
po = [pmat[1,0],pmat[0,1]]
elif result == 10:
po = [pmat[0,1],pmat[1,0]]
else:
po = [pmat[1,1],pmat[1,1]]

return(po)

def strCoop(mode):
if mode == "p1":
return(1)
else:
return(10)

def strDef(mode):
if mode == "p1":
return(0)
else:
return(-10)

def strTfT(lastPlay, mode):
if lastPlay == "D":
if mode == "p1":
return(0)
else:
return(-10)
else:
if mode == "p1":
return(1)
else:
return(10)
Nrep = 10
p1Hist = np.zeros([Nrep,2])
p2Hist = np.zeros([Nrep,2])

for i_rep in range(Nrep):
p1 = strDef("p1")
p2 = strCoop("p2")
po = iterPD(p1,p2,pmat)
p1Hist[i_rep,...] = [p1,po[0]]
p2Hist[i_rep,...] = [p2,po[1]]

lastPlay = "C"
for i_rep in range(Nrep):
p1 = strCoop("p1")
p2 = strTfT(lastPlay, "p2")
if p1 == 0:
lastPlay = "D"
else:
lastPlay = "C"
po = iterPD(p1,p2,pmat)
p1Hist[i_rep,...] = [p1,po[0]]
p2Hist[i_rep,...] = [p2,po[1]]
```

# ケーキの分配ゲーム

```# ケーキの総数
n.cake = 10

# 利得行列
pay.mat = matrix(0,(n.cake+1),(n.cake+1))
for (i.cake in 1:n.cake){
pay.mat[(i.cake+1),1:(n.cake-i.cake+1)] =i.cake
}

# 初期化（各戦略の確率や時間など）
p.cake = runif(n.cake+1)
p.cake = p.cake/sum(p.cake)
max.time = 50
dt = 0.01
t = seq(0,max.time,dt)
n.iter = length(t)
p.hist = matrix(0,nrow = n.iter, ncol = (n.cake+1))
p.hist[1,] = p.cake

# メインのスクリプト
for (i.time in 2:n.iter){
W = colSums(p.cake*t(pay.mat))
W.ave = sum(outer(p.cake,p.cake)*pay.mat)
p.cake = p.cake + p.cake*(W - W.ave)/W.ave * dt
p.hist[i.time,] = p.cake
}

# 結果の可視化
plot(p.hist[,1],ylim=c(0,1),type='l',lwd=2,ylab = 'Proportion',xlab="time")
for (i.strat in 2:(n.cake+1)){
lines(p.hist[,i.strat],col=i.strat,lwd=2)
}
legend("topleft",paste("request = ",0:10),col=1:(n.cake+1),lwd =2,cex=0.75)
```

# dlls

```import numpy as np
import matplotlib.pylab as plt

h = 1e-4 # 0.0001

for idx in range(x.size):
tmp_val = x[idx]
x[idx] = float(tmp_val) + h
fxhP = f(x) # f(x+h)

x[idx] = tmp_val - h
fxhM = f(x) # f(x-h)
grad[idx] = (fxhP - fxhM) / (2*h)

x[idx] = tmp_val # 値を元に戻す

x = init_x
x_history = []

for i in range(step_num):
x_history.append( x.copy() )

return x, np.array(x_history)

def function_2(x):
return (x[0]**2)/20 + x[1]**2

init_x = np.array([-3.0, 4.0])

lr = 0.9
step_num = 100
x, x_history = gradDesc(function_2, init_x, lr=lr, step_num=step_num)

plt.plot( [-5, 5], [0,0], '--b')
plt.plot( [0,0], [-5, 5], '--b')
plt.plot(x_history[:,0], x_history[:,1], '-o')

plt.xlim(-3.5, 3.5)
plt.ylim(-4.5, 4.5)
plt.xlabel("X0")
plt.ylabel("X1")
plt.show()
```

# MCMC

# JAGS and “simple” GIBBS sampling

```# "simple" GIBBS sampling example
# initialization
n.iter = 10000; sigma = 0.1; counter = 0
z = c(6, 2); N = c(8, 7)
a = c(2, 2); b = c(2, 2);
n.par = 2
th.hist = matrix(0, nrow = n.iter*n.par, ncol = n.par)
theta = runif(2)

# function to calc. prob. move
prob.gibbs <- function(theta, proposed, N, z, a, b, idx){
p.th=dbeta(theta[idx], z[idx]+a[idx], N[idx]-z[idx]+b[idx])
p.pro=dbeta(proposed, z[idx]+a[idx], N[idx]-z[idx]+b[idx])
return(p.pro/p.th)
}

# main loop
for (i.iter in 1:n.iter){
for (i.par in 1:n.par){
proposed = theta[i.par] + rnorm(1, mean=0, sd=sigma)
if (proposed > 1) {proposed = 1}
if (proposed < 0) {proposed = 0}
p.move= min(1, prob.gibbs(theta, proposed, N, z, a, b, i.par))
if (runif(1) < p.move){
theta[i.par] = proposed
}
counter = counter + 1
th.hist[counter, ] = theta
}
}
par(mfrow=c(3,1))
HDI.plot(th.hist[,1])
HDI.plot(th.hist[,2])
plot(th.hist[,1],th.hist[,2], type='o',cex=0.1,xlim = c(0,1),ylim=c(0,1))
par(mfrow=c(1,1))

# with JAGS
# model
model {
for (i in 1:Ntotal) {
y[i] ~ dbern(theta[ coin[ i ] ])
}
for (c in 1:Ncoin) {
theta[ c ] ~ dbeta(2, 2)
}
}
# R command
y = c(rep(1,6), rep(0,2), rep(1,2), rep(0,5))
coin = c(rep(1,8), rep(2,7))
dataList = list(Ntotal=length(y),
y = y,
coin = coin,
Ncoin = 2)
parameters = c("theta")
model = jags.model( "~/research/oka/DBDA/ch07/twoVarMCMC.tex",
data = dataList, n.chains = 3, n.adapt = 1000)
update(model, n.iter=500)
CS = coda.samples(model, variable.names = parameters, n.iter = 10000, thin = 5)
mcmcMat<-as.matrix(CS)
par(mfrow = c(3, 1))
plot(mcmcMat[,1], mcmcMat[,2], type='o',ylim = c(0,1), xlim = c(0,1))
HDI.plot(mcmcMat[,1])
HDI.plot(mcmcMat[,2])
```