0

I am currently working on a project where I must analyze data and find a period for the graph. The data contains outliers. I need a function that will make a line of best fit for the function.

I attempted to simply get a sin graph on the plot, but I could not even do that. Can anyone give me a starting hint?

import os
import pyfits as fits
import numpy as np
import pylab
import random
import scipy.optimize
import scipy.signal

from numpy import arange
from matplotlib import pyplot
from scipy.optimize import curve_fit

filename = 'C:\Users\Ken Preiser\Desktop\Space thing\Snapshots\BAT_70m_snapshot_SWIFT_J1647.9-4511B.lc'
namePortion = filename[-39:]

hdulist = fits.open(filename, 'readonly', None, False) #{unpacks file) name, mode, memorymap, savebackup
data = hdulist[1].data

datapoints = 23310

def sinfunc(a, b, c): #I tried graphing a sinfunction, but it did not work...
    return a*np.sin(bx-c) 

time = data.field('TIME')
time = time / 86400.0

timeViewingThreshold = 10
rateViewingThreshold = .01

rate = np.sum(data['RATE'][:,:4], axis=1)

average = np.sum(rate)/23310

error = data.field('ERROR')
error = np.sqrt(np.sum(data['ERROR'][:,:4]**2, axis=1))

print rate.size,(", rate")
print time.size,(", time")

fig = pylab.figure()
ax = fig.add_subplot(111)

ax.set_xlabel('Time')
ax.set_ylabel('Rate')
ax.set_title('Rate vs Time graph: ' + namePortion)

pylab.plot(time, rate, 'o')
pyplot.xlim(min(time) - timeViewingThreshold, max(time) + timeViewingThreshold)
pyplot.ylim(min(rate) - rateViewingThreshold, max(rate) + rateViewingThreshold)
ax.errorbar(time, rate, xerr=0, yerr=error)
pylab.show()

(the outputs) https://i.stack.imgur.com/WYGgj.jpg

Cody Gray - on strike
  • 239,200
  • 50
  • 490
  • 574
user2251518
  • 33
  • 2
  • 5

1 Answers1

0

You're trying to fit points to the model: y = sin(ax + b). Since you're using linear regression, you need a linear model. So one way to do that is compute arcsin for each point and now compute the linear regression. The model is now: arcsin(y) = ax + b. The regression model gives you a and b which is what you're after. You should be able to test this out pretty quickly in excel, then code it up once the nuances are figured out.

Leo Dagum
  • 81
  • 3