在直方图中绘制最佳拟合曲线

Question

我今天花了一整天的时间，想用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势能曲线。不过，我想在这个直方图上画一条曲线，以便得到这条最佳拟合曲线的函数形式。我看到的所有例子都是将曲线画到定义好的函数上，这显然很简单，但在直方图上画曲线似乎是个隐藏的秘密。任何帮助都非常感谢。谢谢。

顺便说一下，这是我完全凭猜测做的尝试：

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()

数据可视化直方图统计分析编程新手数学建模最佳拟合曲线拟合势能曲线

在直方图中绘制最佳拟合曲线

2 个回答

撰写回答