1
$\begingroup$

I am performing a fit on experimental data with zero degree of freedom but with error bar using NonlinearModelFit:

nlmSimple = NonlinearModelFit[p[[All, 1]], a*Exp[-t/k1], {{a, 3000}, {k1, 600}},   t, Weights -> 1/p[[All, 2, 1]]^2,    VarianceEstimatorFunction -> (1 &), MaxIterations -> 200];

where p[[All, 1]]=$\{\{x1,y1\},\{x2,y2\}\}$ is a couple of data point and p[[All, 2, 1]] the associated error on $y_i$ ( which is basically the square-root of the experimental data).

Of course residuals are zero since degree of freedom is zero (=> the function take exactly the value of the data). Consequences : all the standard tools for errors estimations, confidence level etc. of Mathematica are not working.

But I should be able to extract an error since the errors I am using are not equal to zero.

Is there any Mathematica function/option/argument for doing this or I have to build my own tool (and go through the maths... ) ?

Thank you for your help

edit :

The data are describe the radiocative decay of particles. $x_i$ is the number of surviving particle.

So that mean my $x_i$ folllow the binomiale distribution : $\binom{a}{k} p^k(1-p)^{n-k}$ with $p= exp(-y_i/k1)$.

As in practice $x_i>1000$ is quite large and $p \in [ 0.05, 0.9] $ the distribution of $x_i$ convergence toward normal distribution.

$\endgroup$
  • 1
    $\begingroup$ You'll need to be more specific and (maybe) a specific example would help. Your terminology makes it a bit difficult to figure out what can be done. For example, parameters don't have residuals. Is what you are stating that you know what the residual error variance is? In other words your model might be $y_i=a+b x_i + \epsilon_i$ with $\epsilon_i \sim N(0,\sigma^2)$ with $a$ and $b$ unknown but $\sigma^2$ known? $\endgroup$ – JimB Oct 26 '18 at 17:20
  • $\begingroup$ i edit a bit my post to give more information. Yes you are right about the residual except my model is more : $ y_i= a e^{bx_i}+ϵ_i $ $\endgroup$ – Dalnor Oct 26 '18 at 17:33
  • $\begingroup$ @JimB And yes of course I was talking about residual of my data points, my mistake ! $\endgroup$ – Dalnor Oct 26 '18 at 19:16
1
$\begingroup$

Additional edit: I showed how to obtain standard errors for the parameters but maybe what you mean by "estimation of the error" is a confidence interval for the mean prediction. I've added that below. (If you need 95% prediction intervals, then those would be much wider.)

Edit: With the additional information given in the question, it seems clear that you need to perform a Poisson regression.

It appears that you have two counts of particles each over a different length of time. In the olden days one would take the square root of the count and perform a standard linear regression as the square root transformation stabilizes the variance of Poisson counts.

But in the last 20 years or so that's no longer necessary. Software is readily available to explicitly recognize the integral role of the Poisson distribution. Mathematica has GeneralizedLinearModelFit to do so. (It would be even nicer if Mathematica had GeneralizedLinearMixedModelFit but that might have to wait until version 14.)

The model set-up would be the following:

$$y_i \sim Poisson(\lambda_i)$$

$$\log(\lambda_i)=c_0+c_1 t_i$$

Using your notation $\log(\lambda_i)=\log(a)-t_i/k_1$ so GeneralizedLinearModelFit estimates $c_0=\log(a)$ and $c_1=-1/k_1$. To get estimates of the standard errors of the estimators of $a$ and $k_1$ we'll need to use the Delta Method (a.k.a Propagation of Error method).

(* Set parameters *)
a = 3000;
k1 = 600;

(* Generate two data points where y ~ Poisson(a Exp[-t/k1]) *)
t = {1, 100};  (* Lengths of time *)
SeedRandom[12345];
y = Flatten[RandomVariate[PoissonDistribution[a Exp[-#/k1]], 1] & /@ t]
data = Transpose[{t, y}]

(* Perform a Poisson regression *)
(* y ~ Poisson(λ_t) *)
(* log(λ_t)= c0 + c1 t *)
(* where c0 = Log[a] and c1 = -1/k1 *)
(* This is equivalent to y ~ Poisson(a Exp[-t/k1]) *)
glm = GeneralizedLinearModelFit[data, z, z, ExponentialFamily -> "Poisson"];

(* Get parameter estimates which are estimates of Log[a] and 1/k1 *)
{loga, minusRecipk1} = glm["BestFitParameters"];
(* Estimates of a and k1 from Poisson regression *)
ahat = Exp[loga]
(* 2924.6191888146855 *)
k1hat = -1/minusRecipk1
(* 807.5868297019861 *)

Note that we can get the estimates directly without having to run GeneralizeLinearModelFit:

(* Estimate of a *)
N[y[[1]] (y[[1]]/y[[2]])^(1/(t[[2]] - t[[1]]))]
(* 2924.6191888146923 *)

(* Estimate of k1 *)
N[(t[[2]] - t[[1]])/Log[y[[1]]/y[[2]]]]
(* 807.5868297019673 *)

But we need to run GeneralizedLinearModelFit to get the estimate of the parameter covariance matrix. (Although this could be more directly estimated using LogLikelihood.)

Using the Delta Method one can obtain estimates of the standard errors of ahat and k1hat:

covGLM = glm["CovarianceMatrix"];
seahat = Exp[loga] covGLM[[1, 1]]^0.5
(* 54.662928826594296 *)
sek1hat = (1/minusRecipk1^2) covGLM[[2, 2]]^0.5
(* 177.91390393430407 *)

To obtain approximate confidence intervals for the mean prediction we'll also need to use the Delta method.

se[z_, cov_, mle_] := 
  (D[Exp[c0 + c1 z], {{c0, c1}}].cov.D[Exp[c0 + c1 z], {{c0, c1}}] /. 
  {c0 -> mle[[1]], c1 -> mle[[2]]})^0.5

For this particular sample the standard error for a particular time length $t$ is

$$\left(e^{-0.00338036 t} ((0.673417 t-62.4452) t+3086.14)\right)^{0.5}$$

Show the prediction and confidence intervals:

Show[Plot[{glm[z], glm[z] + 1.96 se[z, covGLM, {loga, minusRecipk1}],
  glm[z] - 1.96 se[z, covGLM, {loga, minusRecipk1}]}, {z, t[[1]], t[[2]]},
  PlotStyle -> {Black, {Dotted, Gray}, {Dotted, Gray}}, Frame -> True,
  FrameLabel -> {"t", "Mean count"},
  PlotLabel -> "Predicted mean and\napproximate 95% Conf. Intervals"],
 ListPlot[data, PlotStyle -> {{Red, PointSize[0.03]}}]]

Predicted mean and approximate 95pct confidence intervals

The estimates of the standard errors at $t_1$ and $t_2$ work out to be simply $\sqrt{y_1}$ and $\sqrt{y_2}$, respectively.

$\endgroup$
  • $\begingroup$ Thank you for your help. From your answer I conclude there is no tool for this. I am surprise to have to use random variable to have a result. Now I think the answer is quite simple, error popagation : $ da = \sqrt{ (\frac{da}{dy1})^2 + \frac{da}{dy2})^2 }$ You sould indeed be skeptical. But on one side my systematic errors are way below the statistical error (and well estimate ) so I can do an almost perfect fit with only pure statistical error. One can show that for a define time of data acquisition less data point is more sensitive (smaller error)...so 2 points is the best. $\endgroup$ – Dalnor Oct 26 '18 at 22:07
  • $\begingroup$ I would not conclude that there is no built-in in Mathematica tool to do this but it does seem unlikely. But if one doesn't know the fixed coefficients, it's usually unlikely (and therefore, unbelievable) that one would "know" the variance. So still color me skeptical. $\endgroup$ – JimB Oct 26 '18 at 22:16
  • $\begingroup$ In fact we do know the variance of $x_i$ since it follow binomial distribution $\binom{a}{\exp(-t/k_1)}$ ( we are speaking about radioactive decay). This can be approximate with Poisson law for short $t$ and so variance is $\sqrt{x_i}$. For long $t$ we need indeed an estimation of $a$ and $k_1$ but with one or two recurrence step this is ok. $\endgroup$ – Dalnor Oct 27 '18 at 10:53
  • 1
    $\begingroup$ But $\sqrt{x_i}$ is not the known variance: that is a known "estimate" of the variance. There's a big difference. I think what you have is Poisson regression with two data points. I'll have to delete or at least modify my answer given this new information. $\endgroup$ – JimB Oct 27 '18 at 17:09
  • $\begingroup$ Oh yes you are right. As I am dealing with large number ($x_i > 1000$) I expect the estimator to be very very close to the real variance. So I think in practice the error on the error would not be that large. Also I do not think Poisson regression is ok since my data does not always fulfill the condition for the convergence toward Poisson distribution. However the condition for convergence toward Normal Distribution are fulfill ($x_1>1000$ and $p = exp(-t_i/k1) \in [0.05,0.90]$). I will edit my inital post for adding details $\endgroup$ – Dalnor Oct 27 '18 at 17:52

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.