11
$\begingroup$

I am fitting curves to my data to extract one parameter. However, I am unsure what the certainty of that parameter is and how I would calculate / express its $95$% confidence interval.

Say for a dataset containing data that exponentially decays, I fit a curve to each dataset. Then the information I want to extract is the exponent $b$. I know the values of $t$ and the value of $a$ I am not interested in (thats a variable that comes from the population, not the process Im trying to model).

I use non-linear regression to fit these parameters. However I dont know how to compute the $95$% confidence interval for any method, so broader answers are welcome as well.

$$f= a\cdot e^{-bt}$$ example data and fit

Once I have my value for $b$, how do I calculate its $95$% confidence interval? Thanks in advance!

$\endgroup$
  • $\begingroup$ How do you fit the data? Is your function transformed so as to fit an OLS? $\endgroup$ – johnny Oct 6 '13 at 10:55
  • $\begingroup$ I see from your comments on the answers that you're actually doing nonlinear least squares. You'd have had good answers more quickly if you'd started with that information. I have at least added a relevant tag. $\endgroup$ – Glen_b Oct 7 '13 at 1:23
  • $\begingroup$ @Glen_b Ah Ill be more complete in the future and add it to the question. I did think about however. With some datasets I use absolute L1 distance and other times still I use linear regression. So I was hoping to get a broad answer. $\endgroup$ – Leo Oct 7 '13 at 6:02
  • $\begingroup$ If you want answers for least squares, L1 regression and nonlinear least squares it would be best to be explicit about that. $\endgroup$ – Glen_b Oct 7 '13 at 6:04
10
$\begingroup$

The problem with linearizing and then using linear regression is that the assumption of a Gaussian distribution of residuals is not likely to be true for the transformed data.

It is usually better to use nonlinear regression. Most nonlinear regression programs report the standard error and confidence interval of the best-fit parameters. If yours doesn't, these equations may help.

Each standard error is computed using this equation:

SE(Pi) = sqrt[ (SS/DF) * Cov(i,i) ]

  • Pi : i-th adjustable(non-constant) parameter
  • SS : sum of squared residuals
  • DF : degrees of freedom (the number of data points minus number of parameters fit by regression)
  • Cov(i,i) : i-th diagonal element of covariance matrix
  • sqrt() : square root

And here is the equation to compute the confidence interval for each parameter from the best-fit value, its standard error, and the number of degrees of freedom.

From [BestFit(Pi)- t(95%,DF)*SE(Pi)]  TO  [BestFit(Pi)+
 t(95%,DF)*SE(Pi)] 
  • BestFit(Pi) is the best fit value for the i-th parameter
  • t is the value from the t distribution for 95% confidence for the specified number of DF.
  • DF is degrees of freedom.

    Example with Excel for 95% confidence (so alpha = 0.05) and 23 degrees of freedom: = TINV(0.05,23) DF equals degrees of freedom (the number of data points minus number of parameters fit by regression)

$\endgroup$
  • $\begingroup$ This is exactly what I needed, thank you! I used lsqcurvefit in Matlab, it does not output the confidence interval or standard error. It gives the Lagrange multipliers (?), the residuals and the squared 2-norm of the residuals. Now with that and your answer I can calculate what I need! $\endgroup$ – Leo Oct 6 '13 at 18:31
3
$\begingroup$

If believe an appropriate model for your data is:

$f = ae^{-bt}$

Then you can take a log transform your response data such that an appropriate model is:

$f' = a' -bt$

with $f' = ln(f)$ and $a' = ln(a)$. The transformed data can be fit using simple linear regression and an estimate for the intercept and slope along with standard errors obtained. If the critical t value and standard error are applied to the parameter estimate, a confidence interval for that parameter estimate can be formed. In R:

# Rough simulated data set.
set.seed(1)
a <- 50; b <- 0.2; n <- 25
x <- 1:n
y <- a*(exp(-b * x))
y <- y + rnorm(n, sd=0.25)
y <- ifelse(y>0, y, 0.1)
plot(x,y)

# Linearise:
y2 <- log(y)
plot(x,y2)

# Fit model to transformed data
model <- lm(y2 ~ x)
summary(model)
confint(model)

# Or:
param <- summary(model)$coefficients[, 1]; se <- summary(model)$coefficients[, 2]
param + qt(0.975, 23) * se
param - qt(0.975, 23) * se

If you are using the model for predicting you should be sure to check that the assumptions of SLR have been met - iid $~N(0,\sigma^2)$.

$\endgroup$
  • $\begingroup$ Ah thanks! A very nice and complete answer! This I can use if I do a linearised fit, which I also sometimes do. I hope you dont mind that I accept Harveys answer, as in this case my question was not about the linearised fit. Still a useful answer though! $\endgroup$ – Leo Oct 6 '13 at 18:23

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.