Plotting bestfit curve to histogram

Question

Ive spent all day trying to plot a simple bestfit curve to a histogram using curve_fit and whatever else there is and ive failed miserably. At the end of my tether ive resorted to posting up my lump of code in the hope that someone can at least point me in the right direction towards plotting this curve. Im a complete newbie to programming so please ignore any poor coding. Heres my code:

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

The position histogram i get is a nice inverted lennard-jones potential. however i want to plot a curve to this in order to get the functional form of this bestfit curve. All examples i see plot curves to defined functions which is so obviously simple, but the art of plotting curves to histograms seems to be a hidden secret. Any help is greatly appreciated. Thanks.

BTW, here is my complete guess of an attempt

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

Oh how embarrassing! It was a silly syntax error, heres the offending code, not the ^ instead of **. Terrible. Anyhow, thanks to all for your input and help!

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

Answer 2

As far as I understand, you want to plot the histogram as a curve instead of the bars? Simply use the data binspos, npos returned by the hist function with the plot function. The problem is that there is one more bin edge than data points, so you'll have to compute the centers of the bins:

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

And then just plot the histogram data:

plt.plot(bin_centers, npos)

If you really want to do curve fitting here, you'll probably have to use the bin_centers as input data for the x axis, hopefully something like this will work:

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