Instantly share code, notes, and snippets.

# petigura/gist:1927518 Created Feb 27, 2012

Linear Least Squares
 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