Linear combination fitting

[Caro]

Hello,

I couldn't yet find a solution for the following task and am wondering on whether it is possible to implement it using lmfit and on how to start. I'm quite new to fitting using python so please be kind.

I got a spectrum and want to have a linear combination with component-spectra getting nonnegative weights for those (the nice method to give for constraints led me to lmfit).

If it's not possible, where would be another place to start searching?

Best regards

Caro*

[Jonathan Gutow]

You could probably make it work, but will have to write a function that uses your component spectra (vectors) to generate simulated spectra based on provided weights. However, this actually sounds like a problem that might  benefit from principle component analysis. It depends on the size of the data set you are considering.

Jonathan

[Matt Newville]

Hi Caro,

Yes, you can use lmfit to do such linear combination fitting.  I do this and support it with a GUI form for my domain.  It is sort of overkill to use lmfit when linear least-squares will work, but from the point of view of a non-linear fit, it is pretty fast. I tend to something like this:

###

import numpy as np
import matplotlib.pyplot as plt

import lmfit

sum_to_one = True  # force weights to sum to one

xdat = np.linspace(0, 100, 201)
y0  = np.sin(xdat/9)
y1  = np.cos(xdat/7)
y2  = np.sqrt(xdat/5)/(xdat+30)
ydat = 0.45*y0 + 0.30*y1 + 0.25*y2 + np.random.normal(size=201, scale=0.01)

# build 2d array of components:
ycomps = np.vstack([y0, y1, y2])
ncomps, npts = ycomps.shape

ls_out = np.linalg.lstsq(ycomps.transpose(), ydat, rcond=-1)
weights = ls_out[0]
print("Least Squares Weights: ", weights)

def lincombo_resid(params, xdata, ydata, ycomps):
    ncomps, npts = ycomps.shape
    resid = -ydata
    for i in range(ncomps):
        resid += ycomps[i, :] * params['c%i' % i].value
    return resid

params = lmfit.Parameters()
for i in range(ncomps):
    params.add('c%i' % i, value=weights[i], min=-0.5, max=1.5)

if sum_to_one:
    expr = ['1'] + ['c%i' % i for i in range(ncomps-1)]
    params['c%i' % (ncomps-1)].expr = '-'.join(expr)

expr = ['c%i' % i for i in range(ncomps)]
params.add('total', expr='+'.join(expr))
result = lmfit.minimize(lincombo_resid, params, args=(xdat, ydat, ycomps))

# show results
lmfit.report_fit(result)
plt.plot(xdat, ydat, label='data')
plt.plot(xdat, ydat+result.residual, label='fit')
plt.legend()
plt.show()

#################

which will print out something like

#################

Least Squares Weights:  [0.44905923 0.29889558 0.23803962]
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 6
    # data points      = 201
    # variables        = 2
    chi-square         = 0.01893520
    reduced chi-square = 9.5152e-05
    Akaike info crit   = -1859.27751
    Bayesian info crit = -1852.67090
[[Variables]]
    c0:     0.44901254 +/- 0.00125397 (0.28%) (init = 0.4490592)
    c1:     0.29886831 +/- 0.00126388 (0.42%) (init = 0.2988956)
    c2:     0.25211915 +/- 0.00227465 (0.90%) == '1-c0-c1'
    total:  1.00000000 +/- 0.00227465 (0.23%) == 'c0+c1+c2'
[[Correlations]] (unreported correlations are < 0.100)
    C(c0, c1) = 0.632

#################

That said, I would echo Jonathan's recommendation.  PCA is very useful.  It can tell you how many variable components there are in a series of spectra (or arrays), and it can help determine if another spectrum can be constructed with those same components, but it does not tell you the weight of a set of "known" components.

--
You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/61521f87-7121-4636-8bde-be6cd6548ff5n%40googlegroups.com.

--

--Matt Newville <newville at cars.uchicago.edu> 630-327-7411