pgcudahy pgcudahy - 3 months ago 12
Python Question

How to improve numpy curve fitting

I've got some data from a solar panel charger of time versus current that I'm trying to fit a curve to. The dataset (with times converted by matplotlib.dates.date2num) is at http://pastebin.com/4FAMbbCJ. I put the times into a list called Time and the current into Current.

import numpy
import matplotlib.pyplot as plt

fit = numpy.polyfit(Time, Current, 10)
fit_fxn = numpy.poly1d(fit)

plt.figure(1)
plt.subplot(211)
plt.plot_date(Time, Current, '-r', xdate=True, ydate=False)
plt.title("Current flow over time")
plt.ylabel("milliAmps")
plt.xlabel("Time")

plt.subplot(212)
plt.plot_date(Time, fit_fxn(Time), '-r', xdate=True, ydate=False)
plt.title("Current fxn with time")

plt.show()


The scatter plot comes out fine, but no matter how many coefficients I try with polyfit, I still get a generally straight line. EDIT: To my eyes the current rises and falls with a peak around noon as would be expected with solar power, but the curve puts maximum current at the earliest time and then falls on a straight line from there. I'd add an image but don't have enough reputation points. I think the error is much more likely in my implementation than in polyfit, I'm just looking to see where I screwed up.

Anyone have any ideas on how to find a better curve?

Answer

It is more a problem of the independent variable. You are using polynomial of degree 10 but that means that you will have something like 4.5e+58 -> your coefficients will have to be very small and you will lose precision.

Try shifting the time variable by its minimum this way.

import numpy
import matplotlib.pyplot as plt
t = Time-Time.min()

fit = numpy.polyfit(t, Current, 10)
fit_fxn = numpy.poly1d(fit)

plt.figure(1)
plt.subplot(211)
plt.plot_date(Time, Current, '-r', xdate=True, ydate=False)
plt.title("Current flow over time")
plt.ylabel("milliAmps")
plt.xlabel("Time")

plt.subplot(212)
plt.plot_date(Time, fit_fxn(t), '-r', xdate=True, ydate=False)
plt.title("Current fxn with time")

plt.show()
Comments