In [108]:

    

#次元の呪いをMCMCで実感する


    

In [109]:

    

%matplotlib inline


    

In [110]:

    

#MCMCとは
#MCMCでπを推定する
from __future__ import division
import math
import random

n = 100000
count = 0

for i in range(n):
    u, v = random.uniform(0,1), random.uniform(0,1)
    d = sqrt((u - 0.5)**2 + (v - 0.5)**2)
    
    if d  < 0.5:
        count += 1
        
ratio = count/n

print ratio *4


    

    




3.14784

In [111]:

    

#上のMCMCをn次元球面で行って同様にπを推定す

def curse(n):
    
    m = 10000
    count = 0
    
    for i in range(m):
        
        x = []
        for j in range(n):
            x.append(random.uniform(0,1))
            
        y =  sum(l**2 for l in x)
        
        z = math.pow(y, 1/2)
        
        if z < 1:
            count += 1
    
    ratio = count / m
    
    return math.pow(ratio* (math.pow(2, n))* (math.gamma(n/2.0 +1)), 2/n)


    

In [112]:

    

curse(2)


    

    
Out[112]:





3.1116

In [113]:

    

#次元の呪いを数値で実感
for i in range(2, 15):
    print curse(i)


    

    




3.1384
3.13559468744
3.13177266097
3.09607696262
3.11219614459
3.2305096123
3.14878347402
3.20008164983
3.09201653526
3.0627831821
2.57997446378
3.09086775634
0.0

In [117]:

    

#グラフで実感
import matplotlib.pyplot as plt

def curse2(n):
    
    m = 10000
    count = 0
    W = []
    
    for i in range(1, m+1):
        
        x = []
        for j in range(n):
            x.append(random.uniform(0,1))
            
        y =  sum(l**2 for l in x)
        z = math.pow(y, 1/2)
        
        if z < 1:
            count += 1
            
        p_ratio = count / i
        
        w =  math.pow(p_ratio* (math.pow(2, n))* (math.gamma(n/2.0 +1)), 2/n)
         
        W.append(w)
    
    plt.plot(W)


    

In [119]:

    

#次元の呪いを実感
#最初の一つは気にしないでください
for i in range(1, 16):
    fig = plt.figure(figsize=(30, 30)) #（縦、横のpixel）
    plt.subplot(5,3,i)
    curse2(i)


    

In [ ]: