Fitting an exponential decay using a convolution integral -

Question

I'm fitting the following data where t: time (s), G: counts, f: impulse function:

 t      G      f 
-7200   4.7     0
-6300   5.17    0
-5400   4.93    0
-4500   4.38    0
-3600   4.47    0
-2700   4.4     0
-1800   3.36    0
 -900   3.68    0
    0   4.58    0
  900   11.73   11
 1800   18.23   8.25
 2700   19.33   3
 3600   19.04   0.5
 4500   17.21   0
 5400   12.98   0
 6300   11.59   0
 7200   9.26    0
 8100   7.66    0
 9000   6.59    0
 9900   5.68    0
10800   5.1     0

Using the following convolution integral:

And more specifically:

Where: lambda_1 = 0.000431062 and lambda_2 = 0.000580525.

The code used to perform that fitting is:

#Extract data into numpy arrays
t=df['t'].as_matrix()
g=df['G'].as_matrix()
f=df['f'].as_matrix()
#Definition of the function
def convol(x,A,B,C):
    dx=x[1]-x[0]
    return A*np.convolve(f, np.exp(-lambda_1*x))[:len(x)]*dx+B*np.convolve(f, np.exp(-lambda_2*x))[:len(x)]*dx+C

#Determination of fit parameters A,B,C
popt, pcov = curve_fit(convol, t, g)
A,B,C= popt
perr = np.sqrt(np.diag(pcov))

#Plot fit
fit = convol(t,A,B,C)
plt.plot(t, fit)
plt.scatter(t, g,s=50, color='black')
plt.show()

The problem is that my fit parameters, A, and B are too low and have no physical meaning. I think my problem is related to the step width dx. It should tend to 0 in order to approximate my sum (np.convolve() corresponds a discrete sum of the convolution product) into an integral.

Answer 1

While this is not an answer, I cannot format code in a comment and so post it here. This code shows how to add bounds to curve_fit. Note that if parameter values are returned at or extremely near the bounds there is likely some other problem.

#Determination of fit parameters A,B,C
lowerBounds = [0.0, 0.0, 0.0] # A, B, C lower bounds
upperBounds = [10.0, 10.0, 10.0] # A, B, C upper bounds
popt, pcov = curve_fit(convol, t, g, bounds=[lowerBounds, upperBounds])

Answer 2

I think the problem is that the convolution calculation is incorrect.

import numpy as np
import scipy.optimize
import matplotlib.pyplot as plt

t = np.array([ -7200, -6300, -5400, -4500, -3600, -2700, -1800, -900, 0, 900, 1800, 2700, 3600, 4500, 5400, 6300, 7200, 8100, 9000, 9900, 10800])
g = np.array([ 4.7, 5.17, 4.93, 4.38, 4.47, 4.4, 3.36, 3.68, 4.58, 11.73, 18.23, 19.33, 19.04, 17.21, 12.98, 11.59, 9.26, 7.66, 6.59, 5.68, 5.1])
f = np.array([ 0, 0, 0, 0, 0, 0, 0, 0, 0, 11, 8.25, 3, 0.5, 0, 0, 0, 0, 0, 0, 0, 0])

lambda_1 = 0.000431062
lambda_2 = 0.000580525

delta_t = 900

# Define the exponential parts of the integrals
x_1 = np.exp(-lambda_1 * t)
x_2 = np.exp(-lambda_2 * t)

# Define the convolution for a given 't' (in this case, using the index of 't')
def convolution(n, x):
    return np.dot(f[:n], x[:n][::-1])

# The integrals do not vary as part of the optimization, so calculate them now
integral_1 = delta_t * np.array([convolution(i, x_1) for i in range(len(t))])
integral_2 = delta_t * np.array([convolution(i, x_2) for i in range(len(t))])


#Definition of the function
def convol(n,A,B,C):
    return A * integral_1[n] + B * integral_2[n] + C

#Determination of fit parameters A,B,C
popt, pcov = scipy.optimize.curve_fit(convol, range(len(t)), g)
A,B,C= popt
perr = np.sqrt(np.diag(pcov))

# Print out the coefficients determined by the optimization
print(A, B, C)

#Plot fit
fit = convol(range(len(t)),A,B,C)
plt.plot(t, fit)
plt.scatter(t, g,s=50, color='black')
plt.show()

The values that I get for the coefficients are:

A = 7.9742184468342304e-05
B = -1.0441976351760864e-05
C = 5.1089841502260178

I don't know if the negative value for B is reasonable or not so I have left it as-is. If you want coefficients that are positive, you can constrain them as shown by James.