绘制fitfit曲线到直方图

Question

我花了整整一天的时间尝试绘制一个简单的bestfit曲线使用curve_fit直方图和其他任何事情都是存在的，而且我很失败。在我的系列文章的最后，我采取了张贴代码的方式，希望有人能够至少指出我绘制这条曲线的正确方向。我是一个完全新手编程，所以请忽略任何不良编码。我的继承人代码：

# Velocity Verlet integrator 

def Verlet(x, V, dt, A): 

    x_new = x + V*dt + (A(x,V,R)*(dt**2)/2) 
    V_new =((2)/2+dt)*(V + ((dt/2)*(((48/x_new**13)-(24/x_new**7)) + ((g*2)**(0.5))*R + ((48/x**13)-(24/x**7)) - g*V + ((g*2)**(0.5))*R))) 
    return (x_new, V_new) 


# Start main program 

# Import required libraries 
import numpy as np 
from numpy import array, zeros 
import random 


mu, sigma = 0, 1 # mean and variance 
S = np.random.normal(mu, sigma, 1000) # Gaussian distribution 

g=10 #gamma damping factor 


# Define the acceleration function 

def A(x,V,R): 

    Acc = (((48/x**13)-(24/x**7)) - g*V + ((g*2)**(0.5))*(R)) 

    return Acc 

# Set starting values for position and velocity 
x = array([3]) 
V = array([0]) 




N = 100000 # number of steps 
M = 10 # save position every M_th step 
dt = 0.01 # interval 

# Lists for storing the position and velocity 
Xlist = zeros([1,N/M]) #define vector dimensions 
Vlist = zeros([1,N/M]) 

# Put the initial values into the lists 
Xlist[:,0] = x 
Vlist[:,0] = V 

# Run simulation 

print "Total number of steps:", N 
print "Saving position and velocity every %d steps." % (M) 
print "Start." 
for i in range(N/M): 
    # Run for M steps before saving values 
    for j in range(M): 
      # Update position and velocity based on the current ones 
      # and the acceleration function 
      R = random.choice(S) # selects a different random number from S each time 
      x, V = Verlet(x, V, dt, A) #calculates new pos and vel from Verlet algorithm 

    # Save values into lists 
    Xlist[:, i] = x 
    Vlist[:, i] = V 
print ("Stop.") 




#Define the range of steps for plot 
L = [] 

k=0 
while k < (N/M): 
    L.append(k) 
    k = k + 1 


#convert lists into vectors 
Xvec = (np.asarray(Xlist)).T #Transpose vectors for plotting purpose 
Vvec = (np.asarray(Vlist)).T 


KB=1.3806488*(10**(-23)) #Boltzmann constant 

AvVel = (sum(Vvec)/len(Vvec)) #computes average velocity 
print "The average velocity is: ", AvVel 


Temp = ((AvVel**2)/(3*KB)) #Temperature calculated from equipartition theorem 
print "The temperature of the system is: ", Temp, 


################# 
##### plots ##### 
################# 


# Plot results 
from matplotlib import pyplot as plt 


#plots 1-d positional trjectory vs timestep 
plt.figure(0) 
plt.plot(Xvec,L) 
# Draw x and y axis lines 
plt.axhline(color="black") 
plt.axvline(color="black") 

plt.ylim(0, N/M) 
plt.xlim(0,6) #set axis limits 

plt.show() 




#plots postion distribution histogram 
plt.figure(1) 
num_bins = 1000 
# the histogram of the data 
npos, binspos, patches = plt.hist(Xvec, num_bins, normed=1, facecolor='blue', edgecolor = "none", alpha=0.5) 
PH = plt.hist(Xvec, num_bins, normed=1, facecolor='blue', edgecolor = "none", alpha=0.5) 
plt.xlabel('Position') 
plt.ylabel('Timestep') 
plt.title(r'Position distribution histogram') 

plt.show() 

位置直方图我得到的是一个漂亮的倒兰纳 - 琼斯势。然而，我想绘制一条曲线来获得这个bestfit曲线的功能形式。所有的例子，我看到绘制曲线到明确简单的定义函数，但曲线绘制直方图的艺术似乎是一个隐藏的秘密。任何帮助是极大的赞赏。谢谢。

顺便说一句，这是我的一个尝试完整的猜测

from scipy.optimize import curve_fit 

def func(x,a): 
    return (-a(1/((x^12)-(x^6)))) 


p0 = [1., 0., 1.] 

coeff, var_matrix = curve_fit(func, binspos, npos, p0=p0) 

# Get the fitted curve 
hist_fit = func(binspos, *coeff) 

plt.plot(binspos, hist_fit, label='Fitted data') 
plt.show() 

Answer 1

据我了解，你要绘制的直方图的曲线，而不是吧？只需使用hist函数返回的数据binspos, npos与plot函数即可。问题是，有比数据点多了一个箱子边缘，所以你必须计算箱的中心：

bin_centers = binspos[:-1] + 0.5 * np.diff(binspos) 

然后就是绘制直方图数据：

plt.plot(bin_centers, npos) 

如果你真的想这样做曲线拟合在这里，你可能不得不使用bin_centers作为输入数据为x轴，希望这样的事情会工作：

coeff, var_matrix = curve_fit(func, bin_centers, npos, p0=p0) 

Answer 2

哦多么尴尬！这是一个愚蠢的语法错误，继承人违规的代码，而不是**而不是**。可怕。无论如何，感谢大家的意见和帮助！

def func(x,a): 
    return (-a(1/((x^12)-(x^6))))