Fitting The Theoretical Equation To My Data

Question

I am very, very new to python, so please bear with me, and pardon my naivety. I am using Spyder Python 2.7 on my Windows laptop. As the title suggests, I have some data, a theoretical equation, and I am attempting to fit my data, with what I believe is the Chi-squared fit. The theoretical equation I am using is

import math

import numpy as np

import scipy.optimize as optimize

import matplotlib.pylab as plt

import csv

#with open('1.csv', 'r') as datafile:
 #   datareader = csv.reader(datafile)
 #   for row in datareader:
  #      print ', '.join(row)

t_y_data = np.loadtxt('exerciseball.csv', dtype=float, delimiter=',', usecols=(1,4), skiprows = 1)


print(t_y_data)

t = t_y_data[:,0]

y = t_y_data[:,1]

gamma0 = [.1]

sigma = [(0.345366)/2]*(len(t))

#len(sigma)

#print(sigma)

#print(len(sigma))

#sigma is the error in our measurements, which is the radius of the object


# Dragfunction is the theoretical equation of the position as a function of time when the thing falling experiences a drag force
# This is the function we are trying to fit to our data
# t is the independent variable time, m is the mass, and D is the Diameter

#Gamma is the value of which python will vary, until chi-squared is a minimum



def Dragfunction(x, gamma):
    print x
    g = 9.8
    D = 0.345366
    m = 0.715
#    num = math.sqrt(gamma)*D*g*x
#    den = math.sqrt(m*g)
#    frac = num/den
#    print "frac", frac

    return ((m)/(gamma*D**2))*math.log(math.cosh(math.sqrt(gamma/m*g)*D*g*t))


optimize.curve_fit(Dragfunction, t, y, gamma0, sigma)

This is the error message I am getting:

return ((m)/(gamma*D**2))*math.log(math.cosh(math.sqrt(gamma/m*g)*D*g*t))
TypeError: only length-1 arrays can be converted to Python scalars

My professor and I have spent about three or four hours trying to fix this. He helped me work out a lot of the problems, but this we can't seem to resolve.

Could someone please help? If there is any other information you need, please let me know.

Answer 1

Your error message comes from the fact that those math functions only accept a scalar, so to call functions on an array, use the numpy versions:

In [82]: a = np.array([1,2,3])

In [83]: np.sqrt(a)
Out[83]: array([ 1.        ,  1.41421356,  1.73205081])

In [84]: math.sqrt(a)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
----> 1 math.sqrt(a)

TypeError: only length-1 arrays can be converted to Python scalars

---

In the process, I happened to spot a mathematical error in your code. Your equation at top says that g is in the bottom of the square root inside the log(cosh()), but you've got it on the top because a/b*c == a*c/b in python, not a/(b*c)

log(cosh(sqrt(gamma/m*g)*D*g*t))

should instead be any one of these:

log(cosh(sqrt(gamma/m/g)*D*g*t))
log(cosh(sqrt(gamma/(m*g))*D*g*t))
log(cosh(sqrt(gamma*g/m)*D*t))     # the simplest, by canceling with the g from outside sqrt

A second error is that in your function definition, you have the parameter named x which you never use, but instead you're using t which at this point is a global variable (from your data), so you won't see an error. You won't see an effect using curve_fit since it will pass your t data to the function anyway, but if you tried to call the Dragfunction on a different data set, it would still give you the results from the t values. Probably you meant this:

def Dragfunction(t, gamma):
    print t
    ...
    return ... D*g*t ...

---

A couple other notes as unsolicited advice, since you said you were new to python:

You can load and "unpack" the t and y variables at once with:

t, y = np.loadtxt('exerciseball.csv', dtype=float, delimiter=',', usecols=(1,4), skiprows = 1, unpack=True)

If your error is constant, then sigma has no effect on curve_fit, as it only affects the relative weighting for the fit, so you really don't need it at all.

Below is my version of your code, with all of the above changes in place.

import numpy as np
from scipy import optimize         # simplified syntax
import matplotlib.pyplot as plt    # pylab != pyplot

# `unpack` lets you split the columns immediately:
t, y = np.loadtxt('exerciseball.csv', dtype=float, delimiter=',',
                  usecols=(1, 4), skiprows=1, unpack=True)

gamma0 = .1 # does not need to be a list

def Dragfunction(x, gamma):
    g = 9.8
    D = 0.345366
    m = 0.715
    gammaD_m = gamma*D*D/m # combination is used twice, only calculate once for (small) speedup
    return np.log(np.cosh(np.sqrt(gammaD_m*g)*t)) / gammaD_m

gamma_best, gamma_var = optimize.curve_fit(Dragfunction, t, y, gamma0)