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)

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)

import numpy as np
import matplotlib.pylab as plt

def numGrad(f, x):
    h = 1e-4 # 0.0001
    grad = np.zeros_like(x)
    
    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 # 値を元に戻す
        
    return grad


def gradDesc(f, init_x, lr=0.01, step_num=100):
    x = init_x
    x_history = []

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

        grad = numGrad(f, x)
        x -= lr * grad

    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()

# 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])