public
Created

Linear Least Squares

  • 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
def PLfit(p0,blin,x,data,model,engine=optimize.fmin_powell,err=1):
"""
Partial Linear Model Fitting.
 
Parameters
----------
p0 : Initial parameters
blin : Boolean array specifying which model parameters are fixed
"""
 
pNL0 = p0[~blin]
 
def obj(pNL):
p = np.zeros(p0.size)
p[~blin] = pNL
p[blin] = llsqfit(pNL,blin,x,data,model)
resid = (data - model(p,x)) / err
cost = (resid**2).sum()
return cost
pNLb = engine(obj,pNL0)
pb = np.zeros(p0.size)
pb[~blin] = pNLb
pb[blin] = llsqfit(pNLb,blin,x,data,model)
return pb
 
def llsqfit(pfixed,blin,x,data,model):
"""
Linear Least Squares Fit
 
Fit model(p,x) to data by linear least squares.
 
Parameters
----------
 
pfixed : These parameters are fixed during the linear fitting.
This is the first argument so that this function can be used with
the scipy.optimize suite of optimizers
blin : vector of length p
[False,True] = [NL param, linear param]
x : Independent variable
data : Data we are trying to fit.
model : Callable. must have the following signature model(p,x).
 
Returns
-------
p1 : Best fit plin from linear least squares
 
"""
nlin = blin.size - pfixed.size
 
# Parameter matrix:
pmat = np.zeros(blin.size)
pmat[~blin] = pfixed
pmat = np.tile(pmat,(nlin,1))
pmat[:,blin] = np.eye(nlin)
 
# Construct design matrix
DS = [model(p[0],x) for p in np.vsplit(pmat,nlin)]
DS = np.vstack(DS)
DS = DS.T
plin = np.linalg.lstsq(DS,data)[0]
return plin

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.