Exponential decay function

Question

I have a model where I assumed a variable X would have an exponential decay over time t. I modeled it in R as

$$N = N_1e^{kt},$$ where $k < 0$ and $t$ is 1:95.

However, looking are some real data that I have captured, it seems that the decay is slightly slower to begin with but is much faster towards the end. The greatest gap in my projected data and the actual data is in the middle.

Beginner question, but I was wondering is that type of decay was a known function, or is there some way for me to model it with the exponential decay function.

I have a plot of the two functions, but am not allowed to upload it here. The plot shows a symmetry between the two sets of data, as if the actual data plot is a mirror image of the model data.

I generate the model data using this

gf <- 1.35 * exp(-0.00946 * 1:95)

and the actual data gf2 is

    c(1.35, 1.33776990615813, 1.32552158334356, 1.31325480792751, 
1.30096935107516, 1.28866497856582, 1.27634145060489, 1.26399852162705, 
1.25163594009013, 1.23925344825924, 1.22685078198044, 1.21442767044339, 
1.20198383593233, 1.18951899356456, 1.17703285101575, 1.1645251082312, 
1.15199545712208, 1.13944358124587, 1.12686915546963, 1.11427184561532, 
1.10165130808571, 1.08900718946957, 1.07633912612483, 1.06364674373805, 
1.05092965685848, 1.03818746840501, 1.02541976914382, 1.01262613713477, 
0.999806137143936, 0.986959320019858, 0.974085222030623, 0.961183364158636, 
0.948253251349736, 0.935294371712896, 0.922306195666412, 0.909288175026071, 
0.896239742030314, 0.883160308296893, 0.870049263704956, 0.856905975195814, 
0.843729785484912, 0.830520011676707, 0.817275943773186, 0.803996843065683, 
0.790681940398468, 0.777330434291129, 0.763941488905233, 0.75051423183888, 
0.737047751730707, 0.723541095652442, 0.70999326626635, 0.696403218720661, 
0.682769857252295, 0.669092031461845, 0.655368532220595, 0.64159808716336, 
0.627779355713794, 0.613910923580398, 0.599991296651462, 0.586018894205227, 
0.571992041337234, 0.557908960489589, 0.543767761945991, 0.529566433130946, 
0.51530282652052, 0.500974645933667, 0.486579430925755, 0.472114538946752, 
0.457577124852233, 0.44296411726136, 0.428272191135991, 0.413497735800733, 
0.398636817423197, 0.383685134710617, 0.368637966230041, 0.353490107290977, 
0.338235793692888, 0.322868608762547, 0.30738136887844, 0.291765980931759, 
0.276013262639391, 0.260112712872577, 0.24405221348003, 0.227817635241112, 
0.211392306411718, 0.19475627881785, 0.177885285875707, 0.160749213513732, 
0.143309764509418, 0.125516708879546, 0.10730146998924, 0.0885651903719612, 
0.0691537518211836, 0.048795218179725, 0.0268838502712772)

Answer 1

First, in its current form, this isn't an exponential. Here is the code with some plots:

gf2 <- c(1.35, 1.33776990615813, 1.32552158334356, 1.31325480792751, 
1.30096935107516, 1.28866497856582, 1.27634145060489, 1.26399852162705, 
1.25163594009013, 1.23925344825924, 1.22685078198044, 1.21442767044339,
1.20198383593233, 1.18951899356456, 1.17703285101575, 1.1645251082312,
1.15199545712208, 1.13944358124587, 1.12686915546963, 1.11427184561532,
1.10165130808571, 1.08900718946957, 1.07633912612483, 1.06364674373805,
1.05092965685848, 1.03818746840501, 1.02541976914382, 1.01262613713477,
0.999806137143936, 0.986959320019858, 0.974085222030623, 0.961183364158636,
0.948253251349736, 0.935294371712896, 0.922306195666412, 0.909288175026071,
0.896239742030314, 0.883160308296893, 0.870049263704956, 0.856905975195814,
0.843729785484912, 0.830520011676707, 0.817275943773186, 0.803996843065683,
0.790681940398468, 0.777330434291129, 0.763941488905233, 0.75051423183888,
0.737047751730707, 0.723541095652442, 0.70999326626635, 0.696403218720661,
0.682769857252295, 0.669092031461845, 0.655368532220595, 0.64159808716336,
0.627779355713794, 0.613910923580398, 0.599991296651462, 0.586018894205227,
0.571992041337234, 0.557908960489589, 0.543767761945991, 0.529566433130946,
0.51530282652052, 0.500974645933667, 0.486579430925755, 0.472114538946752,
0.457577124852233, 0.44296411726136, 0.428272191135991, 0.413497735800733,
0.398636817423197, 0.383685134710617, 0.368637966230041, 0.353490107290977,
0.338235793692888, 0.322868608762547, 0.30738136887844, 0.291765980931759,
0.276013262639391, 0.260112712872577, 0.24405221348003, 0.227817635241112,
0.211392306411718, 0.19475627881785, 0.177885285875707, 0.160749213513732,
0.143309764509418, 0.125516708879546, 0.10730146998924, 0.0885651903719612,
0.0691537518211836, 0.048795218179725, 0.0268838502712772)

plot(gf2) #Plot gf2
plot(diff(gf2)) #Plot the change in gf2

Notice that the last change in gf2 (see the above plot for diff(gf2)) on the far right hand side is the biggest change. If you fit gf2 directly (probably with a ratio of polynomials), the biggest residual will probably be in the last data point. So, you might try a fit for the difference in gf2 and then sum it up.

Looking at the second plot for diff(gf2), it looks like a logarithmic plot that was simply flipped. Here's some code and a plot to show that:

gf2_rev <- rev(gf2) #Reverse gf2
plot(diff(gf2_rev)) #Plot it

If you fit that third graph, you'll get a simple equation. From there, you'll have to flip it back to get the second plot, and then accumulate the sum to get the first plot.

The alternative is to fit a ratio of polynomials that will have a hard time with that far right data point.

Edit 1 (01/14/2012) ================================================

After some playing around with the data, I came up with the following.

First, as described above, that far right data point makes this a little harder, so I decided to fit the difference in gf2 (diff(gf2)) and then sum the results up to get gf2. Here are the results (the fit is in red):

ind <- 1:94
#Here's the fit equation    
diff_gf2_trial <- -0.1283/log(3.666e4 + (2.01*(ind^2)) - (575.4*ind))      

plot(diff(gf2), main="Fit the difference so the far right data point works")
lines(diff_gf2_trial, col="red") #Overlay the fit equation in red

gf2_trial <- cumsum(c(gf2[1], diff_gf2_trial)) #Build the cumulative sum
plot(gf2, main="Sum up the fit to get the estimated gf2")
lines(gf2_trial, col="red") #Overlay the cumulative sum in red

Hope that helps.

Edit 2 (01/14/2012) =========================================

The general form of the equation came from plotting and rough fitting the data in various ways. My target was a clean fit with a model that was as simple as possible. I settled on a negative inverted exponential. The curve looked reasonably quadratic or cubic. After some quick fits, it turned out that a quadratic was good enough. Here's the code for that plot (I realize that the above description sounds like a bunch of weasel-words, however it comes from experience and a lot of playing around with the data).

plot(exp(-1/diff(gf2)), main="Plot of exp(ScaleValue/diff(gf2))")
lines(exp(-0.99/diff(gf2)), col="green")
lines(exp(-0.98/diff(gf2)), col="orange")
lines(exp(-0.97/diff(gf2)), col="purple")

At this point, all I have to do is figure out the ScaleValue for exp(ScaleValue/diff(gf2)) and then perform a linear fit for the quadratic terms that fit the above plot. Here's the code:

library(DEoptim) #An optimizer

#Build a function that allows the Scale Value (scaval) to vary and
#return 1 - Adjusted R-Squared value (so it can be minimized).
#In other words, for each "scaval" value that is tried, fit all of
#the polynomial terms.

minfun <- function(scaval) {
  1 - summary(lm(exp(as.numeric(scaval)/diff(gf2)) ~ poly(ind, degree=2, raw=TRUE)))$adj.r.squared
}

ind <- 1:length(diff(gf2)) #Make an index for the difference in gf2 to be used in minfun
junk <- DEoptim(minfun, lower=-1, upper=-0.0001) #Find the Scale
scaval <- as.numeric(junk$optim$bestmem) #Get the answer returned from the optimizer

#Build the resulting model
res <- lm(exp(scaval/diff(gf2)) ~ poly(ind, degree=2, raw=TRUE))
summary(res)

#Here's the fit equation
diff_gf2_new <- scaval/log(res$coefficients[1] + (res$coefficients[2]*ind) + (res$coefficients[3]*(ind^2)))
plot(diff(gf2), main="Fit the difference so the far right data point works") #Plot the difference in gf2
lines(diff_gf2_new, col="blue") #Overlay the fit equation in blue

Notice that the coefficients for diff_gf2_new may be a little different than diff_gf2_trial. That's because I originally used weights to make sure that the far right data point was represented cleanly. It turns out, those weights aren't really needed. Here's the final fit:

summary(res)

Call:
lm(formula = exp(scaval/diff(gf2)) ~ poly(ind, degree = 2, raw = TRUE))

Residuals:
     Min       1Q   Median       3Q      Max 
-117.090  -64.603   -6.481   57.591  292.487 

Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)                         2.671e+05  2.443e+01   10935   <2e-16 ***
poly(ind, degree = 2, raw = TRUE)1 -4.865e+03  1.187e+00   -4099   <2e-16 ***
poly(ind, degree = 2, raw = TRUE)2  2.161e+01  1.211e-02    1786   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 77.27 on 91 degrees of freedom
Multiple R-squared:     1,  Adjusted R-squared:     1 
F-statistic: 4.742e+07 on 2 and 91 DF,  p-value: < 2.2e-16