2017-09-08 4 views
1

このコードは確率分布psi_0_x_squaredを生成します。 次に、この確率に従ってマルコフ連鎖シミュレーションを実行します。この確率psi_0_x_squaredは、実際には、エネルギーレベルn = 0の位置xにある確率である。この確率に従ってxを1000回動かした後、私は位置xのヒストグラムを生成したい。 (位置周波数)このコードを変更して、選択したデフォルトの「クール」ではなく、高さに応じた色のヒストグラムを出力することができます

''' Markov-chain Monte Carlo algorithm for a particle in a Gaussian potential, 
using the Metropolis algorithm. ''' 
import math, matplotlib.pyplot as plt, random 

def probability(x): 

    #wavefunction n=0 evaluated at position x 
    psi_0_x=math.exp(-x ** 2/2.0)/math.pi ** 0.25 

    #probability n=0 to be at position x 
    psi_0_x_squared= psi_0_x**2 

    return psi_0_x_squared 

data_x=[0] 

x = 0.0  #starts at position 0 
delta = 0.5 #stepsize 

for k in range(1000): #for this number of trials 
    x_new = x + random.uniform(-delta, delta) #I displace it a distance delta 


    if random.uniform(0.0, 1.0) < probability(x_new)/probability(x): 
     x = x_new 
    data_x.append(x) 

#histogram 
cm = plt.cm.get_cmap('cool') 
n, bins, patches= plt.hist(data_x, bins=100, normed=True, color='k') 
bin_centers = 0.5 * (bins[:-1] + bins[1:]) 
col = bin_centers - min(bin_centers) 
col /= max(col) 
for c, p in zip(col, patches): 
    plt.setp(p, 'facecolor', cm(c)) 

plt.show() 

答えて

2

変数nにはバーの高さが含まれています。

a = np.random.normal(size=(1000,)) 

cm = plt.cm.get_cmap('cool') 
n, bins, patches= plt.hist(a, bins=100, normed=True, color='k') 
for c, p in zip(n, patches): 
    plt.setp(p, 'facecolor', cm(c)) 

enter image description here

:とても似

for height, p in zip(n, patches): 
    plt.setp(p, 'facecolor', cm(height)) 

:したがって、これはトリックを行う必要があります

関連する問題