这里以荧光漂白恢复曲线计算恢复时间为例介绍python做曲线拟合。

这是一个典型的荧光漂白恢复曲线¹。其中给了一个时间参数，这个参数就是通过曲线拟合来得到的。

Images were processed with the Fiji open-source software, and recovery curves were analysed in Mathematica (Wolfram Research) by fitting them with a single-exponential decay function.

如文中 Methods 部分所述，使用的是 「single-exponential decay function」 来进行拟合的。

除了指数拟合，曲线拟合常规的有多项式拟合，高斯拟合。这些都对应了不同的统计模型。当需要对曲线进行拟合时，就需要搞清楚应该使用什么样的模型来拟合。

这里推荐的做法首先是多看文献，在本领域内大家对于这类表征数据都使用什么模型来拟合，这个一般是属于比较固定和被广泛接收的方法。所以你看得多了就知道什么表征的数据使用什么方法拟合。但这里还只是「知其然不知其所以然」。

统计学基础

如果希望更进一步，就需要了解这背后的统计分布（事实上属于本科生课程），例如下面几种最常见的：

• 二项分布（Binomial distribution）：是N次伯努利实验成功次数的离散概率分布。
• 泊松分布（Poisson distribution）：是二项分布中阳性率特别小的一种情形， 泊松分布适合于描述单位时间内随机事件发生的次数。参数λ是单位时间(或单位面积)内随机事件的平均发生次数。
• 指数分布（Exponential distribution）：常用于各种寿命的分布。指数分布没有「记忆」，例如一个电子元件正常使用100小时的概率，和它已经使用了200小时后，至少还能再用100小时的概率是相同的。
• 正态分布（Normal distribution）：最常见的分布，数据呈现钟形曲线，大部分数据聚集在平均值附近。

所以这里对于荧光漂白恢复，适用指数分布，但需要进行一些变形：

$ I = 1 - e^{-\lambda t} $

当 t = 0 的时候，归一化强度为 0， 然后当t为无穷大的时候，强度趋于1，这样才跟曲线的情况吻合。

半衰期（half-life，记为T1/2）的概念是指从t时刻开始后的一段时间Δt，期间衰变了原有的一半，如：T1/2=(ln 2)/λ≈1.4427/λ≈0.693/λ

直接从图中还原数据点坐标

将截图保存，然后使用 ImageJ 打开，建立好 x 和 y 方向的像素单位到数值的映射关系，以及确定坐标轴原点如下：

• x: 230>10
• y: 130>0.4
• O: (208, 584)

然后使用 MultiPoint 选择工具点选曲线上的点，然后添加到 roiManager，点击 measure，获得所有数据点的像素坐标，并且保存为 csv 表格文件。然后再使用 python 进行转换。

import pandas as pd
data = pd.read_csv("points.csv")
data['x'] = (data['X']-208)/230*10 - 5
data['y'] = (584 - data['Y'])/130*0.4 + 0.2
# 注意图像的y方向，像素坐标从上往下增大，所以这里是减去
data2 = data.filter(items=['x', 'y'])
data2.to_csv("locs.csv", index=None)

读取复原坐标并进行曲线拟合

将复原的坐标临时保存为 locs.csv 文件之后，在此基础上可进行曲线拟合，看看我们计算出来的 half-time 是否和文献报道的一致。

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

data = pd.read_csv('locs.csv')

plt.figure(figsize=(6,2))
sns.lineplot(x='x', y='y', data=data)
plt.show()

可以看到复原的坐标作图时候和文章中大致相同：

注意这个曲线，其实之后 0 之后的一段才是我们要拟合的数据部分。所以还需要对数据进行切片。由于前面的数据量比较小，可以先 data.head() 方法看看在哪个位置切割。最后选择的是 x>-0.5 的数据。

data2 = data[data['x']>-0.5]

然后，由于文献中并未明确提供拟合函数，我只是根据曲线变化趋势，对「single-exponential decay function」 进行了变形。这个可以试试拟合的效果如何。

from scipy.optimize import curve_fit
import numpy as np

def func(x, a):
    y = 1 - np.exp(-a*x)
    return y

popt, pcov = curve_fit(func, data2['x'], data2['y'])
# popt中为拟合的参数，pcov是参数对应的协方差

y2 = func(data2.x, *popt)
# 代入拟合好的参数到原函数并计算拟合数值
plt.plot(data2.x, data2.y)
plt.plot(data2.x, y2)
plt.show()

可以看到，这个函数的拟合效果并不好：

优化拟合并计算拟合误差

蓝色为原数据，橙色为拟合曲线。这里对拟合函数做进一步调整：

def func2(x, a, b, c):
    y = b*(1 - np.exp(-a*x))+c
    return y

我这里增加了一个缩放因子 b，还有一个平移量 c，使得拟合函数能够更加贴近当前数据的情况。然后再次拟合：

在 jupyter notebook 中打印的参数（按函数输入顺序对应 a，b，c），误差和 R方分别如下：

• 0.20209293 0.50868954 0.39449563
• 0.00688486 0.00909046 0.00919506
• R_squared: 0.9888621337117831

提高拟合效果，除了修改拟合函数，还可以在 curve_fit 函数中增加初始猜测值，拟合方法以及拟合迭代次数等，具体可查阅文档。

所以指数分布的 λ 就是 0.202，按照求半衰期的公式带入，就是 0.693/0.202 = 3.5 秒。这个比图中标注的 1.9 秒要大一些，但在一个数量级内，可能是我复原数据点的时候造成的误差。

利用fsolve求方程近似解

这也有可能是我对于文中的 T1/2 的理解有问题。由于文中给的是相对强度，如果是恢复到 50% 的强度，那么就是要求方程的解。这个时候可以使用 fsolve 的方法，快速计算出一个近似解。

from scipy.optimize import fsolve

def func3(x):
    return func2(x, *popt) - 0.5

x = fsolve(func3, 2)
print(x)

这里构造了一个新函数 func3，当前面的拟合函数 func2 为 0.5 的时候，func3 为0，然后再用 fsolve 求近似解，并且提供了一个猜测的初值。

fsolve 返回的是一个 array，求出来的近似解是 1.15 秒，这个又比文章中的 1.9 秒小。

数据分享

通过百度网盘分享的文件：FRAP.zip