4
$\begingroup$

If I have a data that I fit with NonlinearModelfit that fits a data based on two fitting parameters, c1 and c2.

When I used nlm["ParameterTable"] // Quiet I get the following table:

Image

If I have an equation such as:

eq = (2.303*((70 + 273.15)^2)*(c1/c2))/1000

Is there any code (as opposed to doing it manually) I can use to calculate the value of eq with the combined standard deviation based on the standard deviations of c1 and c2 from the table?.

To clarify I would like to get something like: eq = (2.303*((70 + 273.15)^2)*(8.08318/21.1577))/1000=103.604 but also the standard deviation based on the errors of c1 and c2 as to get something like 103.604 +- standard error

Thank you!

EDIT:

For reference eq comes from:

eq = ((log10q - Log10[qref]) == c1*(Tfp - Tfpref)/(c2 + (Tfp - Tfpref)));
model = Tfp /. Solve[eqn, Tfp][[1]]// FullSimplify;
const = {Tfpref -> 70, qref -> 10/60};
model2 = model /. (const // Rationalize) // FullSimplify;

nlm = NonlinearModelFit[data, {model2, c1 > 5, c2 > 5}, {c1, c2}, 
  log10q];

where everything in eq is known except the fitting parameters c1 and c2

| improve this question | |
$\endgroup$
  • $\begingroup$ Perhaps eq /.nlm["BestFitParameters"]? $\endgroup$ – Ulrich Neumann Jun 25 at 6:17
  • $\begingroup$ Are you wanting eq or an estimate of the standard error for eq? $\endgroup$ – JimB Jun 25 at 6:58
  • $\begingroup$ @JimB what I want is for instance to have the value of eq = (2.303*((70 + 273.15)^2)*(8.08318/21.1577))/1000=103.604 but also the standard deviation based on the errors of c1 and c2 as to get something like 103.604 +- standard error $\endgroup$ – John Jun 25 at 17:31
  • $\begingroup$ @UlrichNeumann using eq /.nlm["BestFitParameters"] only gives me the value of eq (e.g. 103.604) but not its standard error. In other words, it does not give me something like 103.604 +- standard error $\endgroup$ – John Jun 25 at 17:51
7
$\begingroup$

If you're wanting an estimate of the standard error for eq, one approach is to use the Delta Method (a.k.a Propagation of Error if you're in the physical sciences).

(* Data from first example in `NonlinearModelFit` documentation *)
data = {{0, 1}, {1, 0}, {3, 2}, {5, 4}, {6, 4}, {7, 5}};
nlm = NonlinearModelFit[data, Log[c1 + c2 x^2], {c1, c2}, x]

eq = (2.303*((70 + 273.15)^2)*(c1/c2))/1000
(* (271.183 c1)/c2 *)
eq /. nlm["BestFitParameters"]
(* 286.391 *)

f = D[eq, {{c1, c2}}] /. nlm["BestFitParameters"]
(* {190.126, -200.789} *)
se = (f.nlm["CovarianceMatrix"].f)^0.5
(* 261.115 *)

More work but better if the desired function estimator does not have an approximate normal distribution is to use a bootstrap approach.

Addition:

I should note that the true standard error almost certainly doesn't exist as the ratio of two normals have no finite moments. However, the "estimate of the standard error" can (depending on the values of the distributions of the estimators) provide a reasonable confidence interval for the estimate of the ratio (such as in +/- 1.96 standard errors).

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks JimB. When I use the last line of your code I get this error: FittedModel::constr: The property values {CovarianceMatrix} assume an unconstrained model. The results for these properties may not be valid, particularly if the fitted parameters are near a constraint boundary. Is there something I should modify from the code? $\endgroup$ – John Jun 25 at 17:38
  • $\begingroup$ I am using your code like this: eq = (2.303*((70 + 273.15)^2)*((c1 /. nlm["BestFitParameters"])/(c2 /.nlm["BestFitParameters"])))/1000; eq /. nlm["BestFitParameters"]; f = D[eq, {{c1, c2}}] /. nlm["BestFitParameters"]; se = (f.nlm["CovarianceMatrix"].f)^0.5 $\endgroup$ – John Jun 25 at 17:44
  • $\begingroup$ I think perhaps the mistake is in the line f = D[eq, {{c1, c2}}] /. nlm["BestFitParameters"] as I get {0,0} as an output from this line. $\endgroup$ – John Jun 25 at 18:02
  • $\begingroup$ I've included a dataset so now the code is complete (and works as advertised). My guess is that you must have defined c1 and c2 previous to f = D[eq, {{c1, c2}}] /. nlm["BestFitParameters"]. $\endgroup$ – JimB Jun 25 at 19:05
  • $\begingroup$ thank you very much for your help !! $\endgroup$ – John Jun 25 at 19:56
4
$\begingroup$

Assuming the c1/c2 estimates came from a large sample, the eq has a non-central ratio of Gaussians distribution. You'd have to simulate or use NExpectation here:

dc1 = NormalDistribution[8.08318, 0.692171];
dc2 = NormalDistribution[21.1577, 3.13379];
td = TransformedDistribution[(2.303*((70 + 273.15)^2)*(c1/c2))/1000
     , {c1 \[Distributed] dc1, c2 \[Distributed] dc2}];

(* random experiment *)
tdrvts = RandomVariate[td, 1000000];
Histogram[tdrvts]
StandardDeviation[tdrvts] (* around 19.2 *)

(* attempt a near-exact mean and stddev *)
meanEst = NExpectation[x, x \[Distributed] td] (* 106.046 *)
mseEst = Sqrt[Quiet@NExpectation[(x - meanEst)^2, x \[Distributed] td]] (* 19.4 *)
| improve this answer | |
$\endgroup$
  • $\begingroup$ None of the moments exist for the ratio of two normals. So using NExpectation might end up giving one much too large a value depending on the values of the parameters of the two normal distributions. However, despite that condition the "estimated variance" can many times be used to construct reasonably accurate confidence intervals. $\endgroup$ – JimB Jun 25 at 17:14
  • $\begingroup$ @JimB The 19.4 looks like an overestimate compared to the random experiment, so you're probably right - I've also hushed up the NExpectation with Quiet because it was warning me about convergence. $\endgroup$ – flinty Jun 25 at 17:17
  • $\begingroup$ Warnings are good things. We should listen to them more than what we usually do. (I should heed that advice more often than I do.) $\endgroup$ – JimB Jun 25 at 17:21
  • 1
    $\begingroup$ @JimB indeed "...so the operator pressed the P key to override the warning and proceed anyway..." $\endgroup$ – flinty Jun 25 at 17:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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