Difference between 2D and 1D Confidence Regions #848
Comments
|
As the author of this code I have to say please inverstigate, I do not remember the exact details. For the 2D case, I followed a physics PhD-thesis. I am rather sure the 1D-case is correct, but less so for the 2D-part. |
|
Thanks, @Tillsten! Do you have a reference for the PhD thesis? |
|
@ThiagoReschutzegger @Tillsten Yes, thanks. I think there is pretty good evidence that |
|
We can see that they diverge for a linear model. import matplotlib.pyplot as plt
import numpy as np
from numpy.linalg import inv
from scipy.special import erf
import lmfit
# defition of the problem
x = np.linspace(1, 10, 250)
np.random.seed(0)
p = lmfit.Parameters()
p.add_many(('a', 4.), ('b', 4.))
y = 10 * x + 10 + 0.1 * np.random.randn(x.size)
def residual(p):
return p['a']*x + p['b'] - y
# fit the data
mini = lmfit.Minimizer(residual, p, nan_policy='propagate')
out1 = mini.minimize(method='Nelder')
out = mini.minimize(method='leastsq', params=out1.params)
sigmas = [1, 2, 3]
probs = [erf(s / np.sqrt(2)) for s in sigmas]
# Using conf_interval()
bound_a = np.zeros((len(sigmas), 2))
bound_b = np.zeros((len(sigmas), 2))
for i, sig in enumerate(sigmas):
ci = lmfit.conf_interval(mini, out, sigmas=[sig], trace=False)
bound_a[i] = (ci['a'][0][1], ci['a'][-1][1])
bound_b[i] = (ci['b'][0][1], ci['b'][-1][1])
colors = ['red', 'blue', 'green']
fig, ax = plt.subplots(figsize=(8, 8))
# cross marks the mean
plt.plot(out.params['a'].value, out.params['b'].value, 'x', color='gray')
# contour lines by conf_interval2d()
cx, cy, grid = lmfit.conf_interval2d(mini, out, 'a', 'b', 100, 100)
ax.contour(cx, cy, grid, levels=probs, colors=colors, linestyles='-')
# dashed lines by conf_interval()
for i, c in zip(range(3), colors):
ax.plot([bound_a[i][0], bound_a[i][1], bound_a[i][1], bound_a[i][0], bound_a[i][0]],
[bound_b[i][0], bound_b[i][0], bound_b[i][1], bound_b[i][1], bound_b[i][0]],
color=c, linestyle='--', alpha=0.5)
# axis name
ax.set_xlabel('a')
ax.set_ylabel('b')
plt.show()Here's the result
|
|
After further thinking: Could you check if the projection of the 2D-plot gives you the same bounds? |
|
I suspect the problem is with using OK, if I change out = minimizer.leastsq()
# prob = prob_func(result, out)
prob = out.chisqr # just give chi-square, we'll work out meaning laterand then think of import numpy as np
from lmfit import conf_interval, conf_interval2d
from lmfit.lineshapes import gaussian
from lmfit.models import GaussianModel, LinearModel
from lmfit.printfuncs import report_ci
import matplotlib.pyplot as plt
np.random.seed(12)
x = np.linspace(1, 100, num=501)
y = (gaussian(x, amplitude=83, center=47., sigma=11.)
+ 0.03*x + 4.0 +
np.random.normal(size=len(x), scale=0.5) )
mod = GaussianModel() + LinearModel()
params = mod.make_params(amplitude=100, center=50, sigma=20,
slope=0, intecept=2)
out = mod.fit(y, params, x=x)
ci = conf_interval(out, out, sigmas=[1,2,3])
print(out.fit_report())
print("values from conf_interval():")
report_ci(ci)
nsamples = 100
parx = 'amplitude'
pary = 'sigma' # 'intercept', ....
c_x, c_y, c2mat = conf_interval2d(out, out, parx, pary, nsamples, nsamples)
# make contours of sigma-levels by hand from chisqr and redchi
sigma_mat = 5 * np.ones((nsamples, nsamples))
sigma_mat[np.where(c2mat < out.chisqr + 4**2*out.redchi)] = 4.0
sigma_mat[np.where(c2mat < out.chisqr + 3**2*out.redchi)] = 3.0
sigma_mat[np.where(c2mat < out.chisqr + 2**2*out.redchi)] = 2.0
sigma_mat[np.where(c2mat < out.chisqr + 1**2*out.redchi)] = 1.0
# to put X/Y axes in "(value - best)/stderr" from the simple values from leastsq
scaled_x = (c_x - out.params[parx].value)/out.params[parx].stderr
scaled_y = (c_y - out.params[pary].value)/out.params[pary].stderr
sigma_levels = [1, 2, 3, 4]
colors = ['black', 'red', 'blue', 'green']
fig, ax = plt.subplots(figsize=(8, 8))
ax.contour(scaled_x, scaled_y, sigma_mat,
levels=sigma_levels, colors=colors, linestyles='-')
ax.set_xlabel(f'{parx}: (val - best) / stderr')
ax.set_ylabel(f'{pary}: (val - best)/stderr ')
ax.grid(True, alpha=0.2, color='grey')
plt.show()This prints out reports of Reduced chi-square is within an order of magnitude of 1, but not exactly 1. The contour plot of the "set-by-hand" sigma-levels looks like this:
which shows pretty good agreement of the 4 contour levels compared to the stderr from Doing that for a couple of other variable pairs gives:
and
again, showing that getting the sigma-levels from "increase chisqr by N**2*redchi" seems pretty reliable. I think that is not circular, and I think that the use of the F-test for chi-square is a little weak, but I am not trained as a statistician ("How did I get here?") I don't have a stong opinion on how to fix this: Any thoughts or suggestions? |
|
I also do not know how I ended up here. Hopefully I will have time this weekend to find the references. IRC the f-test is valid for nested models and in the case of conf1d case actually is identical to a chi sq2 distribution, since the difference of the dofs=1. But there is a reason why the probability function is easily changeable and I put used definitions at the top of the documentation... |
Sure! But how would that help? The contours of the 2D plot are calculated using the built-in function from matplotlib, which uses ContourPy's function
Me neither. So everything I say can be wrong somehow. Honestly, I am also not sure how to fix the use of import matplotlib.pyplot as plt
import numpy as np
from numpy.linalg import inv
from scipy.special import erf
import lmfit
def calc_prob_approx(x, mu, V):
s = (x - mu).T @ inv(V) @ (x - mu)
return erf(s / np.sqrt(2))
# defition of the problem
x = np.linspace(1, 10, 250)
np.random.seed(0)
p = lmfit.Parameters()
p.add_many(('a', 4.), ('b', 4.))
# add gaussian noise
y = 10 * x + 10 + np.random.normal(0, 1, x.shape)
def residual(p):
return p['a']*x + p['b'] - y
# fit the data
mini = lmfit.Minimizer(residual, p, nan_policy='propagate')
out1 = mini.minimize(method='Nelder')
out = mini.minimize(method='leastsq', params=out1.params)
sigmas = [1, 2, 3]
probs = [erf(s / np.sqrt(2)) for s in sigmas]
# Using conf_interval()
bound_a = np.zeros((len(sigmas), 2))
bound_b = np.zeros((len(sigmas), 2))
for i, sig in enumerate(sigmas):
ci = lmfit.conf_interval(mini, out, sigmas=[sig], trace=False)
bound_a[i] = (ci['a'][0][1], ci['a'][-1][1])
bound_b[i] = (ci['b'][0][1], ci['b'][-1][1])
# Using the covariance matrix approximation
V = out.covar
mu = np.array([out.params['a'].value, out.params['b'].value])
a_stderr = out.params['a'].stderr
b_stderr = out.params['b'].stderr
f = 2
# create a grid of points
a = np.linspace(mu[0] - f * a_stderr, mu[0] + f * a_stderr, 100)
b = np.linspace(mu[1] - f * b_stderr, mu[1] + f * b_stderr, 100)
A, B = np.meshgrid(a, b)
X = np.array([A, B]).T
# calculate the probability at each point
prob = np.zeros(X.shape[:2])
for i in range(X.shape[0]):
for j in range(X.shape[1]):
prob[i, j] = calc_prob_approx(X[i, j], mu, V)
colors = ['red', 'blue', 'green']
fig, ax = plt.subplots(figsize=(8, 8))
# cross marks the mean
plt.plot(out.params['a'].value, out.params['b'].value, 'x', color='gray')
# contour lines by covariance matrix approximation
ax.contour(A, B, prob, probs, colors=colors, linestyles='-.')
# contour lines by conf_interval2d()
nsamples = 100
cx, cy, grid = lmfit.conf_interval2d(mini, out, 'a', 'b', nsamples, nsamples)
ax.contour(cx, cy, grid, levels=probs, colors=colors, linestyles='-')
# dashed lines by conf_interval()
for i, c in zip(range(3), colors):
ax.plot([bound_a[i][0], bound_a[i][1], bound_a[i][1], bound_a[i][0], bound_a[i][0]],
[bound_b[i][0], bound_b[i][0], bound_b[i][1], bound_b[i][1], bound_b[i][0]],
color=c, linestyle='--', alpha=0.5)
# axis name
ax.set_xlabel('a')
ax.set_ylabel('b')
plt.show()We get this result, which shows what I said. Not sure how can this information be helpful. Sending it here just in case it is.
|
|
Oh what I meant is summing over on variable to marginalize it out. In this case the 1D confidence intervals are fitting rather well.
|
|
@Tillsten Yes, I agree that But I would be willing to either modify @ThiagoReschutzegger I'm not sure I follow the latest code (or maybe "interpret the plots"): could that difference just be the re-scaling of the covariance? As a general statement, I would say that this is all confusing enough that I would definitely recommend printing out numerical values (and including those here) and clearly commenting on what you see, how it differs from what you expect, and what you think is the problem. Thanks. |
|
The example added with #852 might help clarify how a map of chi-square values can be used to show the n-sigma contour plots. |
|
I suspect I might be missing something here, and in that case I apologize, but how exactly is this behavior different from the known property that 1d confidence intervals for a single variable are not equal to the projection of 2d confidence intervals for the same confidence level? I'm thinking of this classical figure below (from Numerical Recipes). Again, it's entirely possible that I misunderstood the original issue, in which case please just ignore this comment.
|
|
@mdaeron I believe that Numerical Recipes section is meant to caution the reader to take into account the correlations between parameters when assigning confidence levels -- we do that (or "try to do that") here. If I understand the figure correctly (without actually having re-read the text!), 68% of all the values (the dots) fall between the two horizontal lines bracketed by "68% confidence interval on a1" label, and similarly for a0. One might be tempted to interpret that range as the 68% confidence level for the value of a1. Doing so ignores the bias (correlation) of a1 with a0, and is universally cautioned against. The 68% of (a0, a1) values with the lowest fit residual (chi-square or some other figure of merit) lie within the ellipse, and the full extent of that curve for a1 (and a0) should be used as the 68% / 1-sigma confidence level. All of the estimates of uncertainties in lmfit (including the simplest and automatically estimated values) do take such correlations into account. The issue here is that the boxes in @ThiagoReschutzegger calculated as the n-Sigma levels by Using the code at #852 to simply return the matrix of chi-square, and modifying @ThiagoReschutzegger example to use # contour lines by conf_interval2d()
nsamples = 100
## cx, cy, grid = lmfit.conf_interval2d(mini, out, 'a', 'b', nsamples, samples)
## ax.contour(cx, cy, grid, levels=probs, colors=colors, linestyles='-')
# here grid is a matrix of "chi2_out - out.chisqr", the change in chi-square away from the best fit.
cx, cy, grid = lmfit.conf_interval2d(mini, out, 'a', 'b', nsamples, samples, chi2_out)
sigma = np.sqrt(abs(grid)/out.redchi) # grid matrix of sigma values
probs = scipy.special.erf(sigma/np.sqrt(2)) # matrix of probabilities (needs scipy.special.erf)
ax.contour(cx, cy, probs, levels=probs, colors=colors, linestyles='-')then the boxes (found with
|
OK, thanks for the explanation. I suspect there is an interesting (and valid) “cultural” difference here, as by contrast I was taught (my background being physics/chemistry/applied maths) that there are many ways of arbitrarily defining X % confidence regions in the parameter space. Some of them are silly, of course, but iso-χ2 (the ellipses shown here) and unidimensional “bands” in the parameter space (the vertical or horizontal bands above) all appear to be potentially useful. Sometimes you're really only interested in estimating one parameter and the iso-χ2 approach is not relevant. As far as I can tell this is consistent with the message in that chapter of NumRecipes, which IIRC views the iso-χ2 approach as natural, but not the only valid option. My views above are also heavily influenced by Inverse Problem Theory and Model Parameter Estimation (A. Tarantola, SIAM, 2005). But I realize that this is not a universal point of view.
Thanks, I learned something new! Because I don't use |
|
@mdaeron Most of the people writing and using lmfit have backgrounds in physics or chemistry or related fields. The iso-chi-square view may not be universal, but it is the only view that is relevant here. Anyone who wants to do something different will need to use other tools that give incorrect results.
Not only does I can certainly believe there are people who want to not do that. Some of those people may write textbooks. I have no idea what the reference you cite says, but if it says that there is ever a justification for ignoring correlations between variables, I am sorry you had to study from that book - you were misinformed. |
|
Hey again, Using the code at #852 applied to the first example, we get the result:
Which is perfect. I'm closing the issue. Thanks all! |
- includes example and documentation update Closes: #848
|
Thank for the awesome work, @newville! |
|
@newville, to address the factual part of your answer (because obligatory xkcd) in a civil manner: It's not correlation, it's dimensionality. Consider the iso-χ2 confidence regions for 1,2,3,4... independent Gaussian variables. This will yield greater and greater 1-D confidence intervals as dimensionality increases. (You may argue that is The Right Way, but if variables are independent, what is the rational basis for considering them jointly? Why not add yet another, unrelated variable?) It's not like there is no literature on the subject. That'll be all from me now, I don't want to turn this into a game of ping-pong. |
The uncertainties and confidence levels for every variable must take into account the potential correlations between all variables. To this "by hand" (which
I am not arguing, there is no argument to be had, it is the Right Way. If some variables are actually independent of one another, there is no problem accounting for correlations. Indeed the process will show what that correlation is. Again, all fits done with lmfit (which is, like a decade old at this point) do calculate uncertainties that account for correlations and report these values -- the fit report sorts correlations and reports those.
I don't know of anyone here who is not aware of that. Indeed the question is why the different ways to measure confident intervals were giving results, not how to account for correlations. This was not changed.
Great to hear that. Cheers. |
Yes, I read the instructions and I am sure this is a GitHub Issue.
Apparently, there is a difference in calculating the confidence region using the
conf_intervalor theconf_interval2dmethod. It was expected that when plotting the regions calculated by both methods (as I did below), the 1D confidence region would be overlapping the maximum points and minima of the 2D region. However, the results show a significant deviation from this expectation.I would like to investigate this issue further, but before doing so, I need to know whether this behavior is expected or not. If it is not expected, I can focus on understanding the underlying reasons behind this behavior.
Can anyone provide any insights or guidance regarding this issue? Any help would be greatly appreciated.
The text was updated successfully, but these errors were encountered: