3
$\begingroup$

I am trying to figure out what FittedModels`ParameterEllipsoid's parameters mean (I want to construct a filled ellipse)

Here is a minimum working example:

data = {{0, 0}, {1, 1}, {2, 2}, {3, 1}, {4, 0}};
model = a + b x^2;
nlfit = NonlinearModelFit[data, model, {a, b}, x];

nlfit["ParameterConfidenceRegion"]

returns:

FittedModels`ParameterEllipsoid[{1.0069,-0.0344828},{2.60148,0.215202},{{-0.996379,0.0850249},{-0.0850249,-0.996379}}]

Graphically spelunking to find what those parameters mean:

Graphics[
 {
  nlfit["ParameterConfidenceRegion"],
  PointSize[0.01],
  Blue,
  Text["P1", p1 = {1.0068, -0.0344}],
  Orange,
  Text["P2", p2 = {2.6014, 0.215202}],
  Red,
  Text["P3", p3 = {-0.9963, 0.0850}],
  Magenta,
  Text["P4", p4 = {-0.0850, -0.9963}]      
  },
 Frame -> True, AspectRatio -> 1]

ellipse and points

The first point is likely the ellipse center. The third point might be one of the foci.

I tried investigating where the second and fourth points come from. They look like normals because their magnitude is about 1–but nothing obvious pops up, and Normalize[(p3 - p1)] doesn't return p4 or p2.

What do the parameters of FittedModels`ParameterEllipsoid[a, b, {c, d}] mean?

$\endgroup$
1
  • 2
    $\begingroup$ I would bet that {{-0.996379, 0.0850249}, {-0.0850249, -0.996379}} is a rotation matrix. $\endgroup$
    – JimB
    Dec 11, 2023 at 17:08

2 Answers 2

6
$\begingroup$

Addition: At the end I've added how the ellipse is created in a more direct manner.

Copying from @Akku14 's answer (How to get rotate angle θ of rotated ellipse?), the following might be what's needed:

data = {{0, 0}, {1, 1}, {2, 2}, {3, 1}, {4, 0}}
model = a + b x^2
nlfit = NonlinearModelFit[data, model, {a, b}, x]
pcr = nlfit["ParameterConfidenceRegion"]

eli[x_, y_, a_, b_] = x^2/a^2 + y^2/b^2 - 1 == 0
rotate = Thread[{x, y} -> pcr[[3]] . {x, y}]
Show[ContourPlot[Evaluate[eli[x, y, pcr[[2, 1]], pcr[[2, 2]]] /. rotate /. 
  Thread[{x, y} -> {x - pcr[[1, 1]], y - pcr[[1, 2]]}]], 
  {x, -2, 4}, {y, -0.4, 0.3}, 
  FrameLabel -> (Style[#, Bold, 18, Italic] &) /@ {"a", "b"}],
ListPlot[{pcr[[1]]}]]

95% confidence ellipse for parameters a and b

Direct calculation of an approximate confidence ellipse for any two parameters in a model

n = Length[data];  (* Number of data points *)
p = 2;  (* Number of parameters in model *)
i1 = 1;  (* Indices of two parameters for joint confidence ellipse *)
i2 = 2;  (* Note in this example there are only two parameters *)
cov = nlfit["CovarianceMatrix"];  (* Covariance matrix for all parameters *)
(* Get the covariance matrix of the two parameters of interest *)
cov = {{cov[[i1, i1]], cov[[i1, i2]]}, {cov[[i1, i2]], cov[[i2, i2]]}};
estimates = nlfit["BestFitParameters"];  (* Parameter estimates *)

(* Ellipse for all confidence levels *)
ellipse = ((({a, b} - {aTrue, bTrue}) . Inverse[cov] . ({a, b} - {aTrue, bTrue})) /. estimates);

(* Contour value for approximate confidence ellipse *)
confLevel = 0.95;  (* 95% confidence level *)
c = InverseCDF[FRatioDistribution[2, n - p - 1], confLevel];

Show[ContourPlot[ellipse, {aTrue, -2, 4}, {bTrue, -0.4, 0.3}, 
  Contours -> {c}, ContourShading -> None, 
  FrameLabel -> (Style[#, Bold, 18, Italic] &) /@ {"a", "b"}],
 ListPlot[{{a, b} /. estimates}]]

95% confidence ellipse

This also works for LinearModelFit and GeneralizedLinearModelFit. Something nearly identical can be done for a confidence ellipsoid for 3 parameters. For 4 parameters you need magic mushrooms or drop acid to visualize the ellipsoid.

$\endgroup$
1
  • $\begingroup$ That’s it. Thanks. $\endgroup$ Dec 11, 2023 at 18:51
7
$\begingroup$

It is indeed (the inverse of) a rotation matrix, as mentioned by @JimB.

{center, semiaxes, matrix} = List @@ nlfit["ParameterConfidenceRegion"];

Graphics[{
  Yellow, GeometricTransformation[Ellipsoid[{0, 0}, semiaxes], {Inverse@matrix, center}], 
  Red, ellipse}, Axes -> True
]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.