Nonlinear fit with LMFIT (parameters values and variance problem)

[Hiram Lucatero]

Hello, I'm new using LMFIT and I'm trying to make a nonlinear fit, but when I get the parameter's values they don't match with the reported values on literature, I've been trying to change the paratemer's range but it seems that python doesn't take them into consideration because only the wrong values appears  with its variance.  

Here is my code:

import lmfit
from lmfit import minimize, Minimizer, conf_interval, conf_interval2d
import numpy as np
import scipy as sp
import matplotlib
import matplotlib.pyplot as plt




z = np.array([0.100, 0.170, 0.179, 0.199, 0.270, 0.320, 0.352, 0.400,
              0.440, 0.480, 0.570, 0.593,
              0.600, 0.680, 0.730, 0.781,
              0.875, 0.880, 0.900, 1.037, 1.300, 1.363, 1.430, 1.530,
              1.750, 1.965, 2.340, 2.360])
H = np.array([69.0, 83.0, 75.0, 75.0, 77.0, 79.2, 83.0, 95.0, 82.6, 97.0
              , 100.3, 104.0, 87.9, 92.0, 97.3, 105.0, 125.0, 90.0, 117.0,
              154.0, 168.0, 160.0, 177.0, 140.0, 202.0, 186.5, 222.0,
              226.0])


Hin = np.array([12.0, 8.0, 4.0, 5.0, 14.0, 5.6, 14.0, 27.0, 7.8, 62.0,
                3.7, 13.0, 6.1, 8.0, 7.0, 12.0, 17.0, 40.0, 23.0,
                20.0, 17.0, 33.6, 18.0, 14.0, 40.0, 50.4, 7.0, 16.0])




p = lmfit.Parameters()
p.add_many(('Ho', 67., True, None, None),
           ('Om', 0.28, True, 0, 0.3),
           ('Ol', 0.7, True, 0, 1))


def fitCurv(p):
    model= p['Ho']*np.sqrt(p['Om']*(1 + z)**3 + p['Ol']
                           + (1 - p['Om'] - p['Ol'])*(1 + z)**2)
    resids= model - H
    weighted = sp.sqrt(resids ** 2 / Hin ** 2)
    return weighted


mini = lmfit.Minimizer(fitCurv, p)
result = mini.minimize()


print(lmfit.fit_report(result.params))

I get these wrong values (they should be Om=0.28 and Ol=0.72)

[[Variables]]
    Ho:   66.8695791 +/- 3.727759 (5.57%) (init= 67)
    Om:   0.18820806 +/- 0.056660 (30.11%) (init= 0.28)
    Ol:   0.48262270 +/- 0.215481 (44.65%) (init= 0.7)
[[Correlations]] (unreported correlations are <  0.100)
    C(Ho, Ol)                    =  0.909 
    C(Om, Ol)                    =  0.822 
    C(Ho, Om)                    =  0.554 

if I try to make a small range, like this:

p.add_many(('Ho', 67., True, None, None),
           ('Om', 0.28, True, 0, 0.3),
           ('Ol', 0.7, True, 0.6, 0.8))

The fit goes weird and doesn't show the variance of the parameters and It seems like  the fit doesn't vary the initial value of "Ol" 0.7. 

[[Variables]]
    Ho:   70.5082253 +/- 0        (0.00%) (init= 67)
    Om:   0.23648318 +/- 0        (0.00%) (init= 0.28)
    Ol:   0.70000000 +/- 0        (0.00%) (init= 0.7)

How can I fix this? The parameters values must be like the reported ones and also I need the parameters variance in order to make the confidence intervals.

Thank you for your time and help.

[Andrew Nelson]

Dear Hiram,

I looked at this, reproduced it, and can't see any problems. Please see https://gist.github.com/andyfaff/be190d9a027512bcd0a6e342f3e11a11. Note that the final step is a Markov Chain Monte Carlo resampling of the problem to get an estimate of the uncertainties. The uncertainties estimated by `emcee` are the same as those obtained you when you used `mini.minimize()`.

I am assuming that `Hin` represents the uncertainties (standard deviation) estimated on `H`.

When you say the "parameter values must be like the reported ones", where were the parameters previously reported?

There are a few ways to go from here:

1) Original analysis was wrong.

2) You have used the wrong equation for the model.

--
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+unsubscribe@googlegroups.com.
To post to this group, send email to lmfi...@googlegroups.com.
Visit this group at https://groups.google.com/group/lmfit-py.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/2bd281a0-0da4-41c2-811f-5635771927d9%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

_____________________________________
Dr. Andrew Nelson


_____________________________________

[Matt Newville]

Hi Hiram,

Just to follow up on Andrews response:

Two questions:

1. If you know the "right answer", why are you taking data and fitting it to a model?  That's only slightly tongue-in-cheek.   I suggest plotting your H(z) and the model with your starting parameters and the model with the refined parameters.  Honestly, I think you'll have a hard time saying one is obviously a better match to your data, especially with the rather large uncertainties.

2. What is `sp.sqrt(resids**2/Hin**2)` supposed to do?   You should return the residual or the weighted residual to be minimized in the least-squares sense.  I don't think it actually matters much for the result here, but why are you taking the norm here?

 

--Matt

[Andrew Nelson]

On 10 May 2017 at 12:54, Matt Newville <newv...@cars.uchicago.edu> wrote:

2. What is `sp.sqrt(resids**2/Hin**2)` supposed to do?   You should return the residual or the weighted residual to be minimized in the least-squares sense.  I don't think it actually matters much for the result here, but why are you taking the norm here?

I was guessing that Hiram just wanted the absolute value, np.abs(resids/Hin) would've sufficed.

[Matt Newville]

Yeah, but that could confuse `leastsq` as it cannot separate negative residuals from positive ones, which it normally expects to see. 

I think for this sort of simple case, it doesn't make an actual difference in the result, but deliberately suppressing the sign of the residual ems like a bad idea in general.

--Matt

[Andrew Nelson]

Good catch. I didn't think about that because I normally use DE. Still I don't think it changes anything here.

--
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+unsubscribe@googlegroups.com.
To post to this group, send email to lmfi...@googlegroups.com.
Visit this group at https://groups.google.com/group/lmfit-py.

To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/CA%2B7ESbrLqdqEqYMfFRrMWAuT7Tc4mYTrgvHqahbyTYV2uTSOzA%40mail.gmail.com.

[Hiram Lucatero]

Hi Matt,

1. If you know the "right answer", why are you taking data and fitting it to a model?  Because I'm learning how to make model fits, and one exercise that would help me with my present work is the one that I describe above, which is data from the Hubble parameter and the model is the LCDM from cosmology. I tried to remake this fitting in order to obtain the values from the matter and energie densities Om and Ol, which 0.72 and 0.28 are the measured values by the telescope Hubble. I started with this because I want to try to get these values but with a more complex model.  

2. What is `sp.sqrt(resids**2/Hin**2)` supposed to do? I searched for information about how to include the data error in the fit and I found this page: http://stackoverflow.com/questions/15255928/how-do-i-include-errors-for-my-data-in-the-lmfit-least-squares-miniimization-an they say that  `sp.sqrt(resids**2/Hin**2)` should give the correctly  weighted fit.

[Hiram Lucatero]

Hi Andrew,

I saw the code  https://gist.github.com/andyfaff/be190d9a027512bcd0a6e342f3e11a11 and it looks pretty good, but I tried to run it on my computer and after 24 hrs python can't obtain the estimate of the uncertainties, maybe I need a better computer... 

The parameters were reported here  https://arxiv.org/pdf/1502.01589.pdf you can find them on page 32,

H0 . . . . . . . . . . . . 67.31 ± 0.96 

ΩΛ . . . . . . . . . . . . 0.685 ± 0.013 

Ωm . . . . . . . . . . . . 0.315 ± 0.013

I was wrong when I said ΩΛ=0.72, Ωm=0.28 hehehe

To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.

To post to this group, send email to lmfi...@googlegroups.com.
Visit this group at https://groups.google.com/group/lmfit-py.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/2bd281a0-0da4-41c2-811f-5635771927d9%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Matt Newville]

Hi Hiram,

On Wed, May 10, 2017 at 1:03 PM, Hiram Lucatero <hiram.l...@gmail.com> wrote:

Hi Matt,

1. If you know the "right answer", why are you taking data and fitting it to a model?  Because I'm learning how to make model fits, and one exercise that would help me with my present work is the one that I describe above, which is data from the Hubble parameter and the model is the LCDM from cosmology. I tried to remake this fitting in order to obtain the values from the matter and energie densities Om and Ol, which 0.72 and 0.28 are the measured values by the telescope Hubble. I started with this because I want to try to get these values but with a more complex model.  

That's all reasonable.  But I also strongly encourage you to plot the data you actually have with your model and the "known values" and your model with the refined values.   When I made such a plot, I could not easily say that one of those two was obviously better than the other.   The fitting algorithm prefers one over the other.  So the question is: are you more prepared to believe the data you have and the model you have, or the "known values"?   Those "known values" probably came from some source, and should have uncertainties associated with them.   Were these "known values" derived from the data you are using?   If not, what difference in the refined parameters would be acceptable?
 

2. What is `sp.sqrt(resids**2/Hin**2)` supposed to do? I searched for information about how to include the data error in the fit and I found this page: http://stackoverflow.com/questions/15255928/how-do-i-include-errors-for-my-data-in-the-lmfit-least-squares-miniimization-an they say that  `sp.sqrt(resids**2/Hin**2)` should give the correctly  weighted fit.

This may come a shock, but sometimes advice on stackoverflow is incorrect.   The proper way to weight a fit by the uncertainties is the data is to return (data-model)/sigma

where sigma is the uncertainty in the data.

Taking the absolute value of the uncertainty is unnecessary as uncertainties cannot be negative.   Taking the absolute value of the residuals suppresses the sign of the residual, which is useful in calculating the partial derivatives.

--Matt