Instantly share code, notes, and snippets.

What would you like to do?
Linear Least Squares
def PLfit(p0,blin,x,data,model,engine=optimize.fmin_powell,err=1):
Partial Linear Model Fitting.
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.
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).
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)
plin = np.linalg.lstsq(DS,data)[0]
return plin
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment