绘制fitfit曲线到直方图
问题描述:
我花了整整一天的时间尝试绘制一个简单的bestfit曲线使用curve_fit直方图和其他任何事情都是存在的,而且我很失败。在我的系列文章的最后,我采取了张贴代码的方式,希望有人能够至少指出我绘制这条曲线的正确方向。我是一个完全新手编程,所以请忽略任何不良编码。我的继承人代码:绘制fitfit曲线到直方图
# 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()
答
据我了解,你要绘制的直方图的曲线,而不是吧?只需使用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)
+0
是的,我可以像你说的那样绘制数据点,但我实际上是在平滑的最佳拟合曲线之后,而不是实际情节本身。我如何指示Python用这些数据绘制这条最佳曲线?感谢您的回复 – Paddy
答
您是否尝试过调用'func(binspos,* p0)'? – silvado
我其实无法做到这一点,无处不在的错误!哇,我真的不认为这么简单的事情可能会如此困难 – Paddy
如果你在你的问题中发布错误信息,人们可能会帮助你。 – silvado