Fitting exponential decays in R, the easy way

September 9, 2018
By

(This article was first published on R on Programming notes, and kindly contributed to R-bloggers)

Exponential decays can describe many physical phenomena: capacitor discharge, temperature of a billet during cooling, kinetics of first order chemical reactions, radioactive decay, and so on. They are very useful functions, but can be tricky to fit in R: you’ll quickly run into a “singular gradient” error. Thankfully, self-starting functions provide an easy and automatic fix. Read on to learn how to use them.

The formula I’ll use in the following examples is:
$$
y(t) \sim y_f + (y_0 – y_f) e^{-\alpha t}
$$

The measured value $y$ starts at $y_0$ and decays towards $y_f$ at a rate $\alpha$. Let’s generate some artificial data so you can replicate the examples:

library(tidyverse)
library(broom)

t = 1:100
y1 = 22 + (53 - 22) * exp(-0.02 * t) %>% jitter(10)
y2 = 24 + (60 - 22) * exp(-0.01 * t) %>% jitter(10)


df <- data.frame(t = t, y = y1, sensor = 'sensor1') %>% 
  rbind(. , data.frame(t = t, y = y2, sensor = 'sensor2'))
sensor1 <- df %>% filter(sensor == 'sensor1')

Our data looks like this:

qplot(t, y, data = df, colour = sensor)

raw_data

Fitting with NLS

nls is the standard R base function to fit non-linear equations. Trying to fit the exponential decay with nls however leads to sadness and disappointment if you pick a bad initial guess for the rate constant ($\alpha$). This code:

nls(y ~ yf + (y0 - yf) * exp(-alpha * t), 
    data = sensor1,
    start = list(y0 = 54, yf = 25, alpha = 1))

fails with:

Error in nls(y ~ yf + (y0 - yf) * exp(-alpha * t), data = sensor1, start = list(y0 = 54,  : 
  singular gradient

Using SSasymp

The solution is to use a self-starting function, a special function for curve fitting that guesses its own start parameters. The asymptotic regression function, SSasymp is equivalent to our exponential decay:

> fit <- nls(y ~ SSasymp(t, yf, y0, log_alpha), data = sensor1)
> fit
Nonlinear regression model
  model: y ~ SSasymp(t, yf, y0, log_alpha)
   data: sensor1
       yf        y0 log_alpha 
   21.884    52.976    -3.921 
 residual sum-of-squares: 0.9205

Number of iterations to convergence: 0 
Achieved convergence tolerance: 8.788e-07

Its formula is a little different from ours, instead of fitting the rate constant $\alpha$ directly:
$$
y(t) \sim y_f + (y_0 – y_f) e^{-\alpha t}
$$

It searches for the logarithm of $\alpha$:

$$
y(t) \sim y_f + (y_0 – y_f) e^{-\exp(\log\alpha) t}
$$

From the fit result, you can plot the fitted curve, or extract whichever other information you need:

qplot(t, y, data = augment(fit)) + geom_line(aes(y = .fitted))

Exponentially decaying data, fitted

For a single curve, it’s easy to guess the approximate fit parameters by looking at the plot, or just trying several values. When fitting many curves however, it is quite convenient to automate the process. Self-starting functions are especially useful combined with dplyr, to fit several experimental conditions in one step.

Here is how we can read out the fit parameters for each sensor in our example data:

df %>% 
  group_by(sensor) %>% 
  do(fit = nls(y ~ SSasymp(t, yf, y0, log_alpha), data = .)) %>% 
  tidy(fit) %>% 
  select(sensor, term, estimate) %>% 
  spread(term, estimate) %>% 
  mutate(alpha = exp(log_alpha))
sensor log_alpha y0 yf alpha
sensor1 -3.920671 52.97573 21.88404 0.019827795
sensor2 -4.611703 61.98803 23.85592 0.009934881

Now we know at one glance the rate constant for each sensor location, or the $y$ value that each position will stabilise at.

For more ideas on how to apply curve fitting with dplyr, check out my previous article on dplyr.

To leave a comment for the author, please follow the link and comment on their blog: R on Programming notes.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Search R-bloggers

Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)