Fit an exponential decay without initial guessing

fit_decay.R

# Function to fit an exponential decay without initial guessing
# Based on https://stackoverflow.com/a/39436209/545346
# Original source:
# Regressions et Equations integrales, Jean Jacquelin
# https://www.scribd.com/doc/14674814/Regressions-et-equations-integrales
# Converted to R by Johan Van de Wauw

fit_decay <- function(x, y)
{
  n= length(x)
  
  y = y[order(x)]
  x = x[order(x)]
  
  dSk = diff(x) * (y[2:n] + y[1:(n-1)]) /2
  Sk = diffinv(dSk)
  
  dx = x - x[1]
  dy = y - y[1]
  
  m1 = matrix(NA,2,2)
  m1[1,1] = sum(dx**2)
  m1[1,2] = m1[2,1] = sum(dx*Sk)
  m1[2,2] = sum(Sk**2)
  
  m2 = c( sum(dx*dy), sum(dy*Sk))
  
  nc = solve(m1) %*% m2
  
  c = nc[2,1]
  
  m3 = matrix(NA,2,2)
  m3[1,1] = n
  m3[1,2] = m3[2,1] = sum(exp(c*x))
  m3[2,2] = sum(exp(2*c*x))
  
  m4 = c(sum(y), sum(y*exp(c*x)))
  
  ab = solve(m3) %*% m4
  
  return (c(ab[1,1], ab[2,1], c))  
}

# test

x = c(0,0.026,0.052,0.079,0.105,0.131,0.157,0.184,0.211,0.237,
      0.263,0.29,0.315,0.342,0.369,0.395,0.4222,0.448,0.473,0.5)
y = c(4.5,4.42,4.11,4.01,3.56,3.32,3.14,3.32,3.22,2.99,2.96,2.54,
      2.46,2.44,2.58,2.42,2.61,2.41,2.19,2.52) 

# test with non-ordered x values
randorder = order(runif(length(x)))
x <- x[randorder]
y <- y[randorder]

decay <- fit_decay(x,y)

plot(x,y)
lx = seq(0,0.5,0.01)
ly = decay[1]+ decay[2] * exp(decay[3]*lx) 
lines(lx,ly)