我花了一整天时间尝试使用 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()
我得到的位置直方图是一个很好的反转 lennard-jones 势。但是我想为此绘制一条曲线以获得该最佳拟合曲线的函数形式。我看到的所有示例都将曲线绘制为定义的函数,这显然很简单,但是将曲线绘制为直方图的艺术似乎是一个隐藏的 secret 。任何帮助是极大的赞赏。谢谢。
顺便说一句,这是我对一次尝试的完整猜测
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()
最佳答案
据我了解,您想将直方图绘制为曲线而不是条形图?只需将 hist
函数返回的数据 binspos, npos
与 plot
函数一起使用即可。问题是 bin 边缘比数据点多一个,因此您必须计算 bin 的中心:
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)
关于python - 绘制直方图的最佳拟合曲线,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/22845396/