「基礎からのベイズ統計学」を読んで突っかかった点3 : 6章

これまで

6章：正規分布モデルにおける基本的な推定

MH法、HMC法のことを考えすぎてベイズ推定の基礎がどこかに行っていたため全然はじめは理解できませんでした。
3章を再度見直したところ理解できました。

f(\theta|x) \propto f(x | \theta) f(\theta)

だけ思い出せれば十分です。

カタログ刷新問題

文中にある重要な仮定は２つあります。

1. 20人の顧客の購入金額$X$は平均$\mu$,分散$\sigma$に従う
2. 事前分布$f(\mu),f(\sigma)$は一様分布と考える。

ここから、尤度関数$f(x|\mu,\sigma)$を求めると、正規分布の積になるので以下になります。

f(x|\mu,\sigma) = \prod_i^n \frac{1}{\sqrt{2 \pi \sigma^2 }}exp(-\frac{(x_i-\mu)^2}{2 \sigma^2} )

事前分布が一様であることから

f(\mu) \propto const

f(\sigma) \propto const

といえます。

以上より、事後分布のカーネルが導出できます。

f(\mu,\sigma|x) \propto \prod_i^n \frac{1}{\sigma}exp(-\frac{(x_i-\mu)^2}{2 \sigma^2} )

となります。

HMC法で導出

位置関数$h(\mu,\sigma)$を以下の様に計算

h(\mu,\sigma) = - \log f(\mu,\sigma|x) = n \log \sigma+ \frac{\sum (x_i - \mu)^2 }{2 \sigma^2}

非常に計算しやすいですね！

これを、$\mu,\sigma$で微分して力に変換すると

F(\mu) =  \frac{\sum (x_i - \mu)}{\sigma^2 }

F(\sigma) = -\frac{n}{\sigma} + \frac{\sum (x_i - \mu)^2}{ \sigma^3}

が得られます。
この運動方程式に従い、マルコフ遷移することで事後分布を計算します。

Pythonで実装

6章のカタログ刷新問題を解いてみました。

上が平均値の分布、下が標準偏差の分布

mu 2916.8161905460106
sigma 847.6892951700972
upper 2500  0.979
upper 3000  0.33345
base  2500  effect upper 0.8  0.1385

が得られます。教科書の表6.2の数値とほぼ一致します。
プログラムではバーンイン期間を設けていないので、少しズレている可能性はあります。
でも結構いい感じです。
積分計算を飛ばして求めることができるのは便利だと実感できました。

import math 
import numpy as np
import random
import matplotlib.pyplot as plt
import scipy.stats as st

sum_x = 58530
sum_x2 = 182172300
n = 20

def F_m(mu,sigma):
    hoge = (sum_x - n*mu)/(sigma**2)
    return hoge

def F_s(mu,sigma):
    hoge = (sum_x2 - 2*mu*sum_x + n*mu*mu)/(sigma**3) - n/sigma
    return hoge    

def h(mu,sigma):
    hoge = n*math.log(sigma) + (sum_x2 - 2*mu*sum_x + n*mu*mu)/(2*sigma**2)
    return hoge

def H(mu,sigma,p1,p2):
    return h(mu,sigma) + p1*p1/2 + p2*p2/2

m = 2000
s = 100

dt = 0.1
L=100

N = round(L/dt)
cul_rep_num = 0

rnd_m = np.empty([0]) #1次元から配列生成:保存用 平均μ mu
rnd_s = np.empty([0]) #1次元から配列生成:保存用　分散σ sigma

theta_N = 2*10**4

for j in range(theta_N):
    p_m = random.gauss(0,1)
    p_s = random.gauss(0,1)
    
    
    m_start = m
    s_start = s
    pm_start = p_m
    ps_start = p_s
    
    for i in range(N):
        pm_05 = p_m + F_m(m,s)*dt*0.5
        ps_05 = p_s + F_s(m,s)*dt*0.5
        
        m_1 = m + pm_05*dt
        s_1 = s + ps_05*dt
                
        pm_1 =  pm_05 + F_m(m_1,s_1)*dt*0.5
        ps_1 =  ps_05 + F_s(m_1,s_1)*dt*0.5
        
        p_m = pm_1
        p_s = ps_1
        m = m_1
        s = s_1
        
    # 受容するか判定
    r = math.exp( H(m_start,s_start,pm_start,ps_start) - H(m,s,p_m,p_s))
    #print(r)
    
    hoge = random.uniform(0,1)
    if hoge < r:
        rnd_m = np.append(rnd_m,m)
        rnd_s = np.append(rnd_s,s)
        cul_rep_num = cul_rep_num + 1
    else:
        m = m_start
        s = s_start 
        p_m = pm_start
        p_s = ps_start
    
    if j%10**3 == 0:
        print(j)
    
    
# ヒストグラム
fig = plt.figure()

ax = fig.add_subplot(2, 1, 1)
ax.hist(rnd_m, bins=100, histtype='barstacked', ec='black',density=True)

#xx = np.linspace(0, 2.5,501)
#ax.plot(xx, st.gamma(11, 0, 1/13.).pdf(xx))

ax = fig.add_subplot(2, 1, 2)
ax.hist(rnd_s, bins=100, histtype='barstacked', ec='black',density=True)

plt.show()  


# 平均値
mu_eap = np.sum(rnd_m)/ len(rnd_m)
sigma_eap = np.sum(rnd_s)/ len(rnd_s)
print("mu",mu_eap)
print("sigma",sigma_eap)

# 平均値が2500以上の割合
hoge = np.count_nonzero(rnd_m >= 2500)
print("upper 2500 ",hoge/len(rnd_m))

# 平均値が3000以上の割合
hoge = np.count_nonzero(rnd_m >= 3000)
print("upper 3000 ",hoge/len(rnd_m))

# 効果量
mu_0 = 2500
count = 0
for i in range(len(rnd_m)):
    hoge = (rnd_m[i] - mu_0) / rnd_s[i]
    
    if hoge > 0.8:
        count = count + 1
        
print("base ",mu_0," effect upper 0.8 ",count/len(rnd_m))
    

#受容率
print("---- accept num -----")
print(cul_rep_num/(theta_N))

工場機器買い替え問題

カタログ刷新問題との違いは１点のみです。

1. 事前分布$f(\mu)$が一様分布ではなく、正規分布$f(\mu) = N(145,0.1)$

であるとわかっていることです。
ただし、分散に関しては何もわかっていないので先程同様一様分布と考えます。

つまり、

f(\mu) = N(145,0.1)

f(\sigma) \propto const

といえます。

先程同様、事後分布のカーネルを求めます。

f(\theta|x) \propto f(x | \theta) f(\theta)

尤度関数$f(x|\mu,\sigma)$を求めると、正規分布の積になるので以下になります。

f(x|\mu,\sigma) = \prod_i^n \frac{1}{\sqrt{2 \pi \sigma^2 }}exp(-\frac{(x_i-\mu)^2}{2 \sigma^2} )

以上より、事後分布のカーネルが導出できます。

f(\mu,\sigma|x) \propto \{ \prod_i^n \frac{1}{\sigma}exp(-\frac{(x_i-\mu)^2}{2 \sigma^2} ) \}

exp(-\frac{(\mu - 145)^2}{2 \times 0.1})

$0.1^2$ではないことに気をつけてください(1敗)。

位置関数$h(\mu,\sigma)$を以下の様に計算

h(\mu,\sigma) = - \log f(\mu,\sigma|x) = n \log \sigma+ \frac{\sum (x_i - \mu)^2 }{2 \sigma^2} 
+\frac{(\mu-145)^2}{2 \times 0.1}

非常に計算しやすいですね！

これを、$\mu,\sigma$で微分して力に変換すると

F(\mu) =  \frac{\sum (x_i - \mu)}{\sigma^2 } - \frac{(\mu-145)}{0.1}

F(\sigma) = -\frac{n}{\sigma} + \frac{\sum (x_i - \mu)^2}{ \sigma^3}

が得られます。$F(\mu)$の正負の符号に注意してください（2敗）。

HMC法を適用してみます。

Pythonで実装

今回は分散が0.1以上の確率を求めるのが目的である。$\sigma$は標準偏差であるため、サンプルした$\sigma$を2乗してやる必要がある。
これを忘れて原因追求に時間がかかった・・・（3敗///4-33　阪神逆転）

結果は、教科書の値に近い結果が得られた。

mu 144.99170500673344
sigma 0.4096875510127533
variance 0.17017241467236396
upper 0.1  0.9901532278322029
upper 0.15  0.654257724189902
---- accept num -----
0.99525

import math 
import numpy as np
import random
import matplotlib.pyplot as plt
import scipy.stats as st

sum_x = 5799.7
sum_x2 = 840919.141
n = 40

def F_m(mu,sigma):
    hoge = (sum_x - n*mu)/(sigma**2)
    hoge = hoge - (mu-145)/(0.1)
    return hoge

def F_s(mu,sigma):
    hoge = (sum_x2 - 2*mu*sum_x + n*mu*mu)/(sigma**3) - n/sigma
    return hoge    

def h(mu,sigma):
    # if sigma > 0.00001:
    #     hoge = n*math.log(sigma) + (sum_x2 - 2*mu*sum_x + n*mu*mu)/(2*sigma**2)
    # else:
    #     hoge = 10e+10
    hoge = n*math.log(sigma) + (sum_x2 - 2*mu*sum_x + n*mu*mu)/(2*sigma**2)
    hoge = hoge + (mu-145)**2 / (2*0.1)
    return hoge

def H(mu,sigma,p1,p2):
    return h(mu,sigma) + p1*p1/2 + p2*p2/2

m = 145
s = 1

dt = 0.01
L=10

N = round(L/dt)
cul_rep_num = 0

rnd_m = np.empty([0]) #1次元から配列生成:保存用 平均μ mu
rnd_s = np.empty([0]) #1次元から配列生成:保存用　分散σ sigma

theta_N = 2*10**4

for j in range(theta_N):
    p_m = random.gauss(0,1)
    p_s = random.gauss(0,1)
    
    
    m_start = m
    s_start = s
    pm_start = p_m
    ps_start = p_s
    
    for i in range(N):
        pm_05 = p_m + F_m(m,s)*dt*0.5
        ps_05 = p_s + F_s(m,s)*dt*0.5
        
        m_1 = m + pm_05*dt
        s_1 = s + ps_05*dt
                
        pm_1 =  pm_05 + F_m(m_1,s_1)*dt*0.5
        ps_1 =  ps_05 + F_s(m_1,s_1)*dt*0.5
        
        p_m = pm_1
        p_s = ps_1
        m = m_1
        s = s_1
        
    # 受容するか判定
    if s > 0:    #分散が負になることは許容できない    
        r = math.exp( H(m_start,s_start,pm_start,ps_start) - H(m,s,p_m,p_s))
        #print(r)
        
        hoge = random.uniform(0,1)
        if hoge < r:
            rnd_m = np.append(rnd_m,m)
            rnd_s = np.append(rnd_s,s)
            cul_rep_num = cul_rep_num + 1
        else:
            m = m_start
            s = s_start 
            p_m = pm_start
            p_s = ps_start
    else:
        m = m_start
        s = s_start 
        p_m = pm_start
        p_s = ps_start
    
    if j%10**3 == 0:
        print(j)
    
#sigmaは標準偏差なので2乗する
rnd_vari = np.square(rnd_s)


# ヒストグラム
fig = plt.figure()

ax = fig.add_subplot(2, 1, 1)
ax.hist(rnd_m, bins=100, histtype='barstacked', ec='black',density=True)

#xx = np.linspace(0, 2.5,501)
#ax.plot(xx, st.gamma(11, 0, 1/13.).pdf(xx))

ax = fig.add_subplot(2, 1, 2)
ax.hist(rnd_vari, bins=100, histtype='barstacked', ec='black',density=True)

plt.show()  


# 平均値
mu_eap = np.sum(rnd_m)/ len(rnd_m)
sigma_eap = np.sum(rnd_s)/ len(rnd_s)
variance_eap = np.sum(rnd_vari)/ len(rnd_vari)

print("mu",mu_eap)
print("sigma",sigma_eap)
print("variance",variance_eap)

# 分散が0.1以上の割合
hoge = np.count_nonzero(rnd_vari >= 0.1)
print("upper 0.1 ",hoge/len(rnd_s))

# 分散が0.15以上の割合
hoge = np.count_nonzero(rnd_vari >= 0.15)
print("upper 0.15 ",hoge/len(rnd_s))

# 効果量
# mu_0 = 2500
# count = 0
# for i in range(len(rnd_m)):
#     hoge = (rnd_m[i] - mu_0) / rnd_s[i]
    
#     if hoge > 0.8:
#         count = count + 1
        
# print("base ",mu_0," effect upper 0.8 ",count/len(rnd_m))
    

#受容率
print("---- accept ratio -----")
print(cul_rep_num/(theta_N))

代表選考ボーダーライン問題

少し難しいです。
前提条件はカタログ刷新問題と同じです。

1. 事前分布は不明であるため、一様分布を使用
2. サンプルデータ$x_i$は正規分布$N(\mu,\sigma)$に従うと考える

尤度関数、事後分布のカーネルが刷新問題同様に得られるので、HMCを使って、$f(\mu|x_{1 \cdots n}),f(\sigma|x_{1 \cdots n})$の確率密度関数に従うサンプルデータが得られます。

得られたサンプルデータが正規分布に従うと考えて（事後分布の関数から推測可）、各サンプルの0.75四分位数を計算。
Kくんは$N(805,100)$に従って、記録を出すと考えられ、先ほど求めた0.75四分位数を超える確率を計算可能。
各サンプルの0.75四分位数に対して、K君が上回る記録を出す確率は累積密度関数から計算できる。
得られたK君が上回る記録を出す確率のサンプル群の期待値・分散・９５％点を出していけば本節の答えが導ける

Python

平均、標準偏差、K君が四分位数を超える確率の確率分布の順に並べました。

mu 786.1929425451208
sigma 16.055175046478478
quartile 3/4 797.0301857014938
p 3/4 0.7704042745873245
x 3/4 0.95  804.294499071472
x 3/4 0.05  790.8521894678336
p 3/4 0.95  92.14613087087847 %
p 3/4 0.05  52.84718942549493 %
---- accept num -----
0.9999

教科書とは若干数値に差異があるような結果

import math 
import numpy as np
import random
import matplotlib.pyplot as plt
import scipy.stats as st

sum_x = 15724
sum_x2 = 12366438
n = 20

def F_m(mu,sigma):
    hoge = (sum_x - n*mu)/(sigma**2)
    return hoge

def F_s(mu,sigma):
    hoge = (sum_x2 - 2*mu*sum_x + n*mu*mu)/(sigma**3) - n/sigma
    return hoge    

def h(mu,sigma):
    hoge = n*math.log(sigma) + (sum_x2 - 2*mu*sum_x + n*mu*mu)/(2*sigma**2)
    return hoge

def H(mu,sigma,p1,p2):
    return h(mu,sigma) + p1*p1/2 + p2*p2/2

m = 786
s = 10

dt = 0.1
L=100

N = round(L/dt)
cul_rep_num = 0

rnd_m = np.empty([0]) #1次元から配列生成:保存用 平均μ mu
rnd_s = np.empty([0]) #1次元から配列生成:保存用　分散σ sigma

theta_N = 2*10**4

for j in range(theta_N):
    p_m = random.gauss(0,1)
    p_s = random.gauss(0,1)
    
    
    m_start = m
    s_start = s
    pm_start = p_m
    ps_start = p_s
    
    for i in range(N):
        pm_05 = p_m + F_m(m,s)*dt*0.5
        ps_05 = p_s + F_s(m,s)*dt*0.5
        
        m_1 = m + pm_05*dt
        s_1 = s + ps_05*dt
                
        pm_1 =  pm_05 + F_m(m_1,s_1)*dt*0.5
        ps_1 =  ps_05 + F_s(m_1,s_1)*dt*0.5
        
        p_m = pm_1
        p_s = ps_1
        m = m_1
        s = s_1
        
    # 受容するか判定
    r = math.exp( H(m_start,s_start,pm_start,ps_start) - H(m,s,p_m,p_s))
    #print(r)
    
    hoge = random.uniform(0,1)
    if hoge < r:
        if j > 10**3:
            rnd_m = np.append(rnd_m,m)
            rnd_s = np.append(rnd_s,s)
        cul_rep_num = cul_rep_num + 1
    else:
        m = m_start
        s = s_start 
        p_m = pm_start
        p_s = ps_start
    
    if j%10**3 == 0:
        print(j)
    
    
# 四分位数導出
rnd_3_4 = rnd_m +  rnd_s * 0.675

# Kの分布N(805,100)に従うXに対して、X>四分位数3/4となる確率
p_34 = 1 - st.norm( loc = 805 , scale = 10 ).cdf(rnd_3_4)    
quartile_p3_eap = np.sum(p_34)/ len(p_34)
print("p 3/4",quartile_p3_eap)   
    
# ヒストグラム
fig = plt.figure(figsize=(5, 10))

ax = fig.add_subplot(3, 1, 1)
ax.hist(rnd_m, bins=100, histtype='barstacked', ec='black',density=True)

#xx = np.linspace(0, 2.5,501)
#ax.plot(xx, st.gamma(11, 0, 1/13.).pdf(xx))

ax = fig.add_subplot(3, 1, 2)
ax.hist(rnd_s, bins=100, histtype='barstacked', ec='black',density=True)

ax = fig.add_subplot(3, 1, 3)
ax.hist(p_34, bins=100, histtype='barstacked', ec='black',density=True)

plt.show()  

# 平均値
mu_eap = np.sum(rnd_m)/ len(rnd_m)
sigma_eap = np.sum(rnd_s)/ len(rnd_s)
quartile_3_eap = np.sum(rnd_3_4)/ len(rnd_3_4)
print("mu",mu_eap)
print("sigma",sigma_eap)
print("quartile 3/4",quartile_3_eap)

# 四分位数の95%値
sort_34 = np.sort(rnd_3_4)
line = 0.95
hoge = sort_34[round(len(sort_34)*line)]
print("x 3/4 "+str(line)+" ", hoge)

# 四分位数の5%値
sort_34 = np.sort(rnd_3_4)
line = 0.05
hoge = sort_34[round(len(sort_34)*line)]
print("x 3/4 "+str(line)+" ", hoge)

# p3/4の95%値
sort_p34 = np.sort(p_34)
line = 0.95
hoge = sort_p34[round(len(sort_p34)*line)]
print("p 3/4 "+str(line)+" ", hoge*100,"%")

# p3/4の5%値
sort_p34 = np.sort(p_34)
line = 0.05
hoge = sort_p34[round(len(sort_p34)*line)]
print("p 3/4 "+str(line)+" ", hoge*100,"%")
   

#受容率
print("---- accept num -----")
print(cul_rep_num/(theta_N))

体内時計問題

1. 事前分布は一様
2. 各群のサンプルデータは正規分布に従う

カタログ刷新問題と同じで各群の尤度関数から各群の平均・分散事後分布が求められる。
HMCで各群にの事後分布に従う$\mu_1,\mu_2$のサンプルを計算し、平均の差の分布関数が求められる。
差が$c$以上となる確率が計算できる。

Python

平均差の分布

運動時平均 30.835789717521354
安静時平均 31.9817793066478
運動時 sd 1.0914644448746031
安静時 sd 1.7440599327296324
平均差 1.1459895891264498
u > 0  0.9906276326874474
u > 1  0.6245787700084247
---- accept num -----
0.9996

まぁ、だいたい一致

# -*- coding: utf-8 -*-
"""
Created on Thu Sep 14 18:13:28 2023
HMC法検証

正規分布タイプ
体内時計問題
2群差
"""

import math 
import numpy as np
import random
import matplotlib.pyplot as plt
import scipy.stats as st

# X 実験群　運動時
sum_x = 616.71
sum_x2 = 19036.1699
n = 20

# Y 対照群　安静時
sum_y = 639.56
sum_y2 = 20501.8104

def F_m(mu,sigma):
    hoge = (sum_x - n*mu)/(sigma**2)
    return hoge

def F_s(mu,sigma):
    hoge = (sum_x2 - 2*mu*sum_x + n*mu*mu)/(sigma**3) - n/sigma
    return hoge    

def F_my(mu,sigma):
    hoge = (sum_y - n*mu)/(sigma**2)
    return hoge

def F_sy(mu,sigma):
    hoge = (sum_y2 - 2*mu*sum_y + n*mu*mu)/(sigma**3) - n/sigma
    return hoge    


def h(mu,sigma):
    hoge = n*math.log(sigma) + (sum_x2 - 2*mu*sum_x + n*mu*mu)/(2*sigma**2)
    return hoge

def H(mu,sigma,p1,p2):
    return h(mu,sigma) + p1*p1/2 + p2*p2/2

def hy(mu,sigma):
    hoge = n*math.log(sigma) + (sum_y2 - 2*mu*sum_y + n*mu*mu)/(2*sigma**2)
    return hoge

def Hy(mu,sigma,p1,p2):
    return hy(mu,sigma) + p1*p1/2 + p2*p2/2


m = 30
s = 1

my = 33
sy = 1

dt = 0.01
L=10

N = round(L/dt)
cul_rep_num = 0

rnd_m = np.empty([0]) #1次元から配列生成:保存用 平均μ mu
rnd_s = np.empty([0]) #1次元から配列生成:保存用　分散σ sigma

rnd_my = np.empty([0]) #1次元から配列生成:保存用 平均μ mu
rnd_sy = np.empty([0]) #1次元から配列生成:保存用　分散σ sigma

rnd_u = np.empty([0]) #1次元から配列生成:保存用　平均差　ux - uy

theta_N = 2*10**4

for j in range(theta_N):
    p_m = random.gauss(0,1)
    p_s = random.gauss(0,1)
    
    p_my = random.gauss(0,1)
    p_sy = random.gauss(0,1)
    
    
    m_start = m
    s_start = s
    pm_start = p_m
    ps_start = p_s
    
    my_start = my
    sy_start = sy
    pmy_start = p_my
    psy_start = p_sy
    
    for i in range(N):
        
        # X
        pm_05 = p_m + F_m(m,s)*dt*0.5
        ps_05 = p_s + F_s(m,s)*dt*0.5
        
        m_1 = m + pm_05*dt
        s_1 = s + ps_05*dt
                
        pm_1 =  pm_05 + F_m(m_1,s_1)*dt*0.5
        ps_1 =  ps_05 + F_s(m_1,s_1)*dt*0.5
        
        p_m = pm_1
        p_s = ps_1
        m = m_1
        s = s_1
        
        #Y
        pmy_05 = p_my + F_my(my,sy)*dt*0.5
        psy_05 = p_sy + F_sy(my,sy)*dt*0.5
        
        my_1 = my + pmy_05*dt
        sy_1 = sy + psy_05*dt
                
        pmy_1 =  pmy_05 + F_my(my_1,sy_1)*dt*0.5
        psy_1 =  psy_05 + F_sy(my_1,sy_1)*dt*0.5
        
        p_my = pmy_1
        p_sy = psy_1
        my = my_1
        sy = sy_1
        
    # 受容するか判定
    r = math.exp( H(m_start,s_start,pm_start,ps_start) - H(m,s,p_m,p_s))
    ry = math.exp( Hy(my_start,sy_start,pmy_start,psy_start) - Hy(my,sy,p_my,p_sy))
    #print(r)
    
    u = my - m
    
    hoge = random.uniform(0,1)
    if hoge < r and  hoge < ry:        
        if j > 10**3:
            rnd_m = np.append(rnd_m,m)
            rnd_s = np.append(rnd_s,s)
            rnd_my = np.append(rnd_my,my)
            rnd_sy = np.append(rnd_sy,sy)
            rnd_u = np.append(rnd_u,u)
            
        cul_rep_num = cul_rep_num + 1
    else:
        m = m_start
        s = s_start 
        p_m = pm_start
        p_s = ps_start
        
        my = my_start
        sy = sy_start 
        p_my = pmy_start
        p_sy = psy_start
    
    
        
    if j%10**3 == 0:
        print(j)
    
    
# ヒストグラム
fig = plt.figure(figsize=(5, 3))

ax = fig.add_subplot(1, 1, 1)
ax.hist(rnd_u, bins=100, histtype='barstacked', ec='black',density=True)



plt.show()  

# 平均値
mu_eap = np.sum(rnd_m)/ len(rnd_m)
myu_eap = np.sum(rnd_my)/ len(rnd_my)
sigma_eap = np.sum(rnd_s)/ len(rnd_s)
sigma_y_eap = np.sum(rnd_sy)/ len(rnd_sy)
u_eap = np.sum(rnd_u)/ len(rnd_u)
print("運動時平均",mu_eap)
print("安静時平均",myu_eap)
print("運動時 sd",sigma_eap)
print("安静時 sd",sigma_y_eap)
print("平均差",u_eap)

# u>0
hoge = np.count_nonzero(rnd_u >= 0)
print("u > 0 ",hoge/len(rnd_u))

# u>1
hoge = np.count_nonzero(rnd_u >= 1)
print("u > 1 ",hoge/len(rnd_u))

#受容率
print("---- accept num -----")
print(cul_rep_num/(theta_N))