Re: Digest for lmfit-py@googlegroups.com - 1 update in 1 topic

[Mark Sprague]

Thanks, Charles. I will give you a call today.

-Mark

On Tue, Jun 7, 2022 at 9:01 AM <lmfi...@googlegroups.com> wrote:

lmfi...@googlegroups.com

Google Groups

Topic digest
View all topics

• Accounting for input data correlations within lmfit - 1 Update

Accounting for input data correlations within lmfit

Denis Leshchev <leshche...@gmail.com>: Jun 06 05:05PM -0400 Tarek, all,   Yeah, thanks for that. I agree that it seems completely possible to add > way to handle covariances that would lead to "negative effective > uncertainties", which could happen mathematically, but have no clear > meaning.   *Diagonal* elements of covariance matrix, square root of which gives us uncertainties (epsilon in lmfit notation) cannot be negative as you correctly point out! However, correlations between data points can be negative. Those will result in *off-diagonal* elements being negative. An important condition that covariance matrix must satisfy is that it has to be *positive semi definite*, which is not the same as not having any negative elements. So a covariance matrix where off-diagonal elements are negative can be a valid one!     Like, I think the idea is to extend uncertainty `epsilon` to `sqrt(epsilon > one thing to have that work in a carefully crafted example and quite > another to support it for others. And we would want to extend "weight" to > "1/epsilon". ;)   The 2020 thread provides lengthy explanations with working examples showing that this is a bad idea. However, you have the right intuition that we should do something *like* square root of the covariance matrix. Square rooting matrices is not as straightforward as scalars though! The good news is that linear algebra does extend the concept of square roots into the matrix realm/ The square roots of matrices can be computed in many ways (wiki page on matrix square root is a good resource to learn about it) and are not unique except for very special cases. This ambiguity does not matter though In the context of optimization/fitting. In fact, Cholesky decomposition I used in my solution is one of the definitions of matrix square root. I used it for convenience since it requires only one line of code.   Maybe the solution is "well known" and we (or to be clear: I) just don't > morons, you just don't care to solve the problem" actually does not > understand the problem at all. We need solutions in Python that are robust > and that we can support.   In the entire script I enclosed only two lines of code take care of the correlations. If an implementation of two lines of code that use standard libraries (such as numpy) is not robust/simple then I don’t know what is. In fact, scipy’s *curve_fit uses exactly the same approach with cholesky decomposition of covariance matrix:* https://github.com/scipy/scipy/blob/4cf21e753cf937d1c6c2d2a0e372fbc1dbbeea81/scipy/optimize/_minpack_py.py#L762 Their implementation is a bit different technically (they really take advantage of the Cholesky factor's triangular structure) and it most likely results in a better computational performance. However, the methodology is the same.   I have further extended my example to show that curve_fit yields exactly the same result as lmfit with cholesky factor: # imports import numpy as np from scipy import linalg, stats from lmfit import Minimizer, Parameters import matplotlib.pyplot as plt   from scipy.optimize import curve_fit   # set seed to 1 np.random.seed(1)   #%% # create true data   def func(params, x): amp = params['amp'] phase = params['phase'] freq = params['freq'] decay = params['decay'] offset = params['offset'] return amp * np.sin(x*freq + phase) * np.exp(-x*x*decay) + offset   def func_array(x, amp, phase, freq, decay, offset): params = {'amp' : amp, 'phase' : phase, 'freq' : freq, 'decay' : decay, 'offset' : offset} return func(params, x)   params = Parameters() params.add('amp', value=5) params.add('phase', value=-0.1) params.add('freq', value=2.0) params.add('decay', value=0.025) params.add('offset', value=10.0)   x = np.linspace(0, 15, 501) data_true = func(params, x) p0_arr = [5, -0.1, 2.0, 0.025, 10.0] #%% # create noisy data with covariance that has off-diagonal elements corelation_width = 1e-1 base = np.exp(-(x/(np.sqrt(2)*corelation_width))**2) A = linalg.toeplitz(base) C = A @ A.T *1e-2 + np.eye(x.size)*1e-3 data = np.random.multivariate_normal(data_true, C)   # pre-compute some matrices to speed up calculations C_inv = np.linalg.inv(C) L = np.linalg.cholesky(C_inv)   #%%   def residual_chol(params, x, data, L): dy = data-func(params, x) return L.T @ dy # multiplication by Cholesky factor is the key component in handling correlated noise   minimizer = Minimizer(residual_chol, params, fcn_args=(x, data, L)) result = minimizer.minimize() # find residuals residual = result.residual   p_arr, p_arr_cov = curve_fit(func_array, x, data, p0=p0_arr, sigma=C) p_arr_err = np.sqrt(np.diag(p_arr_cov))   # figures plt.figure(1) plt.clf()   plt.subplot(311) plt.plot(x, data, 'k.') plt.plot(x, func(result.params, x), 'r-') plt.plot(x, func_array(x, *p_arr), 'b:') plt.xlabel('x') plt.ylabel('data') plt.title('data/fit')   plt.subplot(312) plt.plot(x, data - func(result.params, x), 'k-') plt.xlabel('x') plt.ylabel('$r$') plt.title('raw residuals')   plt.subplot(313) plt.plot(x, result.residual, 'k-') plt.xlabel('x') plt.ylabel('$L^Tr$') plt.title('normalized residuals')   plt.tight_layout()   print('\t\tlmfit\t\t\tcurve_fit') for i, k in enumerate(params.keys()): print(k, 'value ', '\t', result.params[k].value, '\t', p_arr[i]) print(k, 'stderr', '\t', result.params[k].stderr, '\t', p_arr_err[i])   The result prints out the following: amp value 5.112655779597885 5.112665047701475 amp stderr 0.1295239096465614 0.12952611761549376 phase value -0.08243341946021458 -0.08243446730640175 phase stderr 0.02497059590587641 0.024970471047126793 freq value 1.99514253241486 1.9951432017813735 freq stderr 0.009654727425451334 0.009655074417609226 decay value 0.028541209815362276 0.028541424039505475 decay stderr 0.0016738967218192092 0.001674024833945412 offset value 9.99514278795249 9.995142838118364 offset stderr 0.03720335936338103 0.03720337660121335   It is disappointing to see, especially from a scientist in a closely > related field, but there seems no point in engaging with someone with that > attitude.   It is also disappointing to have conversations where other parties (including scientists from related or unrelated fields) completely dismiss proposed solutions, explanations, code, examples and just state that the person “actually does not understand the problem at all” with no attempts to actually look into what they are trying to say. Despite this I still love lmfit and hope to get to the bottom of this one day! ;-)   Cheers, Denis.     Cheers, Denis.

Back to top

You received this digest because you're subscribed to updates for this group. You can change your settings on the group membership page. To unsubscribe from this group and stop receiving emails from it send an email to lmfit-py+u...@googlegroups.com.