使用卷积积分拟合指数衰减- [英] Fitting an exponential decay using a convolution integral -

问题描述

我正在拟合以下数据，其中t：时间（s），G：计数，f：脉冲函数：

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

使用以下卷积积分：

更具体地说：

其中： lambda_1 = 0.000431062 和 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()

问题是我适合参数A和B太低，没有任何物理意义。我认为我的问题与步长 dx 有关。为了趋近我的总和（ np.convolve（）对应于卷积积的离散总和），它应该趋向于0。

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.

推荐答案

我认为问题是卷积计算不正确。

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()

我得到的系数为：

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

我不知道B的负值是否合理，所以我将其保留照原样。如果您希望系数为正，则可以如James所示约束它们。

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.

这篇关于使用卷积积分拟合指数衰减-的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！