Linear Least Squares

gistfile1.py

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