public
Created

Fitting a partially linear model

  • Download Gist
gistfile1.py
Python
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
def fit1T(pNL0,t,f):
"""
Fit Single transit
"""
dpNL0 = np.array([0.2,0.2])
 
objp = lambda p,t,f : obj1Tlin(p,t,f) + (((p-pNL0)/dpNL0)**2).sum()
pNL = optimize.fmin(objp,pNL0,args=(t,f),disp=False)
pL = linfit1T(pNL,t,f)
pFULL = np.hstack( (pNL[0],pL[0],pNL[1],pL[1:]) )
if pFULL[1] < 0:
pFULL[1] = 0
 
return pFULL
 
 
def obj1Tlin(pNL,t,f):
"""
Single Transit Objective Function. For each trial value of epoch
and width, we determine the best fit by linear fitting.
 
Parameters
----------
pNL - The non-linear parameters [epoch,tdur]
"""
pL = linfit1T(pNL,t,f)
pFULL = np.hstack( (pNL[0],pL[0],pNL[1],pL[1:]) )
model = keptoy.P051T(pFULL,t)
 
resid = ((model - f)/1e-4 )
obj = (resid**2).sum()
return obj
 
 
def linfit1T(p,t,f):
"""
Linear fit to 1 Transit.
 
Depth and polynomial cofficents are linear
 
Parameters
----------
p : [epoch,tdur]
t : time
f : flux
 
Returns
-------
p1 : Best fit [df,pleg0,pleg1...] from linear fitting.
"""
epoch = p[0]
tdur = p[1]
tDS = trendDS(t)
ndeg = 3
# Construct lightcurve design matrix
plc = np.hstack(( epoch,1.,tdur,list(np.zeros(ndeg+1)) ))
lcDS = keptoy.P051T(plc,t)
 
DS = np.vstack((lcDS,tDS))
p1 = np.linalg.lstsq(DS.T,f)[0]
return p1
 
 
def trendDS(t):
ndeg=3
# Construct polynomial design matrix
tDS = []
for i in range(ndeg+1):
pleg = np.zeros(ndeg+1)
pleg[i] = 1
tDS.append( keptoy.trend(pleg,t) )
tDS = np.vstack(tDS)
return tDS

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.