Skip to content

ModelResult.eval_uncertainty does not work for complex models #900

Closed
@newville

Description

@newville

Description

ModelResult.eval_uncertainty is not working properly for models that generate complex data. Running the method gives a warning about casting complex to float, and returns a real array of the same length as the model.

See https://groups.google.com/g/lmfit-py/c/tqaLwSBH-Ts/m/HBDP0y4HEAAJ

A Minimal, Complete, and Verifiable example
import numpy as np
from lmfit import Model

def cmplx(f, omega, areal, aimag, off, sigma):
    return (areal*np.cos(f*omega + off) + 1j*aimag*np.sin(f*omega + off))*np.exp(-f/sigma)

f = np.linspace(0, 10, 501)
dat = cmplx(f, 4, 10, 5, 0.2, 4.5) + (0.1+0.2j)*np.random.normal(scale=0.25, size=len(f))


mod = Model(cmplx)
params= mod.make_params(omega=5, areal=5, aimag=5,
                        off={'value': 0.5, 'min': -2, 'max':2},
                        sigma={'value': 3, 'min': 1.e-5, 'max': 1000})

result = mod.fit(dat,  params=params, f=f)
print(result.fit_report())

dfit = result.eval_uncertainty()
print('Data  : ', len(dat), dat.dtype)
print('Fit   : ', len(result.best_fit), result.best_fit.dtype)
print('D_Fit : ', len(dfit), dfit.dtype)
Fit report:
[[Model]]
    Model(cmplx)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 49
    # data points      = 1002
    # variables        = 5
    chi-square         = 1.59641755
    reduced chi-square = 0.00160122
    Akaike info crit   = -6444.87518
    Bayesian info crit = -6420.32641
    R-squared          = (0.9996115675746108+2.8498175850731404e-05j)
[[Variables]]
    omega:  3.99961585 +/- 2.4279e-04 (0.01%) (init = 5)
    areal:  9.99808550 +/- 0.00755744 (0.08%) (init = 5)
    aimag:  4.99883013 +/- 0.00590546 (0.12%) (init = 5)
    off:    0.20069043 +/- 7.0226e-04 (0.35%) (init = 0.5)
    sigma:  4.50042395 +/- 0.00497052 (0.11%) (init = 3)
[[Correlations]] (unreported correlations are < 0.100)
    C(omega, off)   = -0.7328
    C(areal, sigma) = -0.6982
    C(aimag, sigma) = -0.4351
    C(areal, aimag) = +0.3021
/Users/Newville/Codes/lmfit-py/lmfit/model.py:1602: ComplexWarning: Casting complex values to real discards the imaginary part
  fjac[key][i] = (res1[key] - res2[key]) / (2*dval)
Data  :  501 complex128
Fit   :  501 complex128
D_Fit :  501 float64

Version information

Python: 3.10.10 | packaged by conda-forge | (main, Mar 24 2023, 20:17:34) [Clang 14.0.6 ]

lmfit: 1.2.1.post49+g78c9994, scipy: 1.11.0, numpy: 1.25.0,asteval: 0.9.30, uncertainties: 3.1.7

Link(s)

See https://groups.google.com/g/lmfit-py/c/tqaLwSBH-Ts/m/HBDP0y4HEAAJ

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions