曲线拟合中常数的估计方法

Question

我将得到两个长度相同的数组。 An example is here. 在示例中，时间数组用于 x 轴，强度数组用于 y 轴。我被要求用稍微改变的对数正态概率密度函数拟合图形（域是 $t>0$）：

（我真的不知道如何在Whosebug中使用LaTeX）

\begin{equation*}
f(t) = \frac{a}{\sigma t\sqrt{2\pi}}\exp(\frac{-(\ln(t)-m)^2}{2\sigma^2}) + b
\end{equation*}

换句话说，如果两个数组的图形如下所示：

我想确定拟合常数 $a$、$b$、$\sigma$ 和 $m$，使得 $f(t)$ 类似于：

我的计划是使用scipy.optimize.curve_fit。更准确地说：

intensity = np.load('intensity.npy')
time = np.load('time.npy')

def func(t, a, b, m, s): #s=sigma and n = m, just to recopy my formula for f
    t = np.array(t)
    return b + a*np.exp((-1)*(np.log((t))-n)**2/(2*s**2)) / (t)*s*np.sqrt(2*np.pi)

initial = #something to work out
popt, pcov = curve_fit(func, time, intensity,initial)
plt.plot(time, func(time, *popt), linewidth=0.6, color = 'black')
plt.plot(time,intensity,linewidth=0.1, color = 'black')
plt.xlabel("Time(Seconds)")
plt.ylabel("Average Intensity")
plt.savefig('result.png', dpi = 1200)

如果未给出初始估计，该方法将无法正常工作。然而，这里的问题是阵列强度和时间是给我的。我无法控制初始估计，因此需要一种数学方法来获得对拟合常数的一般估计。这是我卡住的地方，因为我真的不熟悉对数正态分布。有什么方法可以帮助我快速估计常数吗？

Answer 1

这确实不是 scipy 或 Python 问题，也许更适合姊妹网站 https://datascience.stackexchange.com/

然而，这是我提出的用于获得良好初始值的解决方案。它确实依赖于对数正态分布的知识，例如 https://en.wikipedia.org/wiki/Log-normal_distribution

您尝试拟合的曲线是经过偏移（b）和缩放（a）的对数正态概率密度。首先我们要估算这些。由于对数正态密度取值在 $[0,+infty)$ 中，参数 b 的初始估计可以设置为 intensity 的最低值。此外，定义对数正态密度的积分应为 1，因此 a 应设置为 intensity 偏移 b 的积分。之后我们在 intensity 中有一个（真实的）对数正态密度，我们可以根据分布的均值及其方差估计 m 和 s （实际上我们只使用第二个时刻） .这个逻辑在下面的函数get_initial_values中实现。我还修改了你的 func 以确保 t 始终严格为正，否则它会弄乱日志，导致被零除等。它还有一个类型 n vs m.这是完整的代码

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
def func(t, a, b, m, s): #s=sigma and n = m, just to recopy my formula for f
    t = np.array(t)
    t = np.maximum(t,1e-4) # guard against zeros or negative values
    return b + a*np.exp((-1)*(np.log(t)-m)**2/(2*s**2)) / (t)*s*np.sqrt(2*np.pi)

def get_initial_values(t,x):
    min_x = np.min(x)
    integral_x = np.sum((x[:-1] - min_x)*np.diff(t))
    
    a = integral_x
    b = min_x
    normalized_x = (x - b)/a

    mean_t = np.sum((normalized_x[:-1] *t[:-1] )*np.diff(t))
    mean_t2 = np.sum((normalized_x[:-1] *t[:-1]**2 )*np.diff(t))
    

    s = np.sqrt(np.log(mean_t2) - 2*np.log(mean_t))
    m = np.log(mean_t) - s*s/2

    return np.array([a,b,m,s])

initial = get_initial_values(time, intensity)
popt, pcov = curve_fit(func, time, intensity,initial)

此处绘制了您的数据结果：

曲线拟合中常数的估计方法

Method for estimating a constant in curve fitting

optimization

curve-fitting

scipy