asymmetric least squares

als.py

from scipy import sparse
from scipy.sparse.linalg import spsolve
import numpy as np

def als(self, y, lam, p, niter):
        '''Asymmetric Least Squares Smoothing algorithm of Eilers and Boelens (2005)

        From Stackoverflow 29185844 by Sparrowcide as modified by jpatina
        
        Good choices for the parameters are:

        - lam: 1000000
        - p: .0001
        - niter: 10

        ref: Eilers, P. H. C. and Boelens, H. F. M. [Baseline Correction with Asymmetric Least Squares Smoothing](https://zanran_storage.s3.amazonaws.com/www.science.uva.nl/ContentPages/443199618.pdf)  Leiden University Medical Centre, October 2005

        Example:
            lam = 1000000
            p = .0001
            niter = 10
            raw = [1,2,3,3,5,6,5,5,4,5,3,3,2,3,2,2,1,1,1]
            baseline = als(raw, lam=lam, p=p, niter=niter)
            print(baseline)

            > [1.00218168 1.00208511 1.00198852 1.00189194 1.00179535 1.00169875
            1.00160214 1.00150552 1.00140889 1.00131225 1.0012156  1.00111894
            1.00102228 1.0009256  1.00082892 1.00073224 1.00063555 1.00053886
            1.00044217]

        Args:
            y: A list of response values from set of inputs in ascending order of time etc.
            lam: lambda value, try varying between 10^2 and 10^9
            p: asymmetry parameter. Try varying between 0.001 and 0.1
        
        Return: list of baseline values
        '''
    L = len(y)
    D = sparse.diags([1,-2,1],[0,-1,-2], shape=(L,L-2))
    w = np.ones(L)
    for i in range(niter):
        W = sparse.spdiags(w, 0, L, L)
        Z = W + lam * D.dot(D.transpose())
        z = spsolve(Z, w*y)
        w = p * (y > z) + (1-p) * (y < z)
    return z