42、随机扰动与相关函数分析

随机扰动与相关函数分析

1. 概率密度函数估计

在研究随机扰动时，概率密度函数（pdf）的估计是一个活跃的研究领域。当 μ = 4 时，二次映射可以生成与随机过程产生的反正弦 pdf 难以区分的分布函数。而且，有研究证明，对于随机过程产生的每个 pdf，都存在无限多个混沌确定性过程可以产生相同的 pdf。

目前，估计 pdf 的方法有多种。一种可能是使用基于人工神经网络或核估计的方法。而当下非常流行的方法是基于特征函数 ϕ( f)。通常，先计算特征函数再进行傅里叶变换，比使用传统方法获得概率密度更容易。同时，估计经验特征函数的方法也在不断发展，并且已经有了一些可用的软件包，如 MATLAB 的 STABLE 软件包和 R 的 gmm 包。

2. 矩的估计

从有限数据集估计矩的效果如何呢？通常测量均值或方差的方法是对同一实验进行多次试验。对于每次试验，计算矩 εi 的值，经过 N 次试验后，计算 ε 的期望值 E[ε]：
[E[\varepsilon] = \frac{1}{N} \sum_{i=1}^{N} \varepsilon_i]
这就是 ε 的样本均值。如果 E[ε] 等于真实值 ˆε（例如通过从 pdf 计算矩得到），则称估计 E[ε] 是无偏的。例如，计算平均值来估计均值就是对一阶矩的无偏估计。反之，如果 E[ε] 不等于 ˆε，则存在偏差 b[ε]：
[b[\varepsilon] = E[\varepsilon] - \hat{\varepsilon}]
当 b[ε] 不等于 0 时，估计就是有偏的。

总估计误差由均方误差给出，即 E[(εi - ˆε)²]。估计可能有偏差，因为均方误差包含非系统性随机分量和系统性偏差分量。就像射击手的子弹集中在靶心右侧，存在明显的偏差，但如果子弹孔之间的间距小，随机分量也会小。

统计学家们致力于开发无偏方法来估计各种矩。例如，当均值为零时，样本值 xi 方差的无偏估计为：
[E[(x - E[x])^2] {unbiased} = \frac{1}{N - 1} \sum {i=1}^{N} x_i^2]
而有偏估计为：
[E[(x - E[x])^2] {biased} = \frac{1}{N} \sum {i=1}^{N} x_i^2]
虽然通常认为无偏估计总是更好，但实际并非总是如此，因为可能无法同时最小化均方误差的两个分量，所以需要进行权衡。近年来对最大似然估计器的重视，就是为了开发能最小化均方误差的方法。

3. 信号平均

黑箱对输入的响应常常被不相关的噪声波动所掩盖。信号功率 Psig 与噪声功率 Pnoise 的比值称为信噪比（SNR）：
[SNR = \frac{P_{sig}}{P_{noise}} = \frac{E[x_{sig}^2]}{E[x_{noise}^2]} = \frac{E[x_{sig}^2]}{\sigma^2}]
这里假设噪声是方差为 σ² 的白噪声。

信号平均技术用于减少记录中噪声的影响，同时增强所需信号的贡献。这些技术将重复的、相同的刺激应用于“开环”动态系统，主要有两个目标：一是测量刺激施加与系统响应之间的时间延迟或潜伏期；二是表征系统对施加刺激的响应。在临床神经生理学中使用诱发电位就是一个研究得很好的例子。

对于 n 个采样信号 z 的白噪声功率，可以很容易地估计为：
[E\left[\frac{1}{n}E\left[\sum_{i=1}^{n} x_i^2\right]\right] = \frac{1}{n^2} E\left[\sum_{i=1}^{n} x_i^2\right] = \frac{1}{n^2} \sum_{i=1}^{n} E_i[x^2] = \frac{1}{n^2} n\sigma^2 = \frac{1}{n}\sigma^2]
由此可见，随着 n 的增加，E[x²noise] 减小。

如果对刺激的响应频率高于其奈奎斯特频率，即过采样。在这种情况下，带宽受限的信号可以完全恢复，并且 SNR 会随着 n 的增加而增加：
[\frac{E[x^2]}{\frac{1}{n}\sigma^2} = \frac{nE[x^2]}{\sigma^2} = n\frac{P_{sig}}{P_{noise}} = nSNR]

为了改进信号平均过程，有以下几种常用技术：
- 时间锁定 ：通常将信号与刺激进行时间锁定，最简单的方法是使用刺激的开始来触发数据收集系统。
- 奇偶试验平均 ：研究人员经常分别对奇数和偶数试验进行平均。奇数和偶数时间平均值的差异可以估计噪声功率，它们的和则是使用所有试验的信号功率。
- 噪声要求 ：为了使信号平均产生预期的结果，噪声必须是不相关的。在实验室和临床环境中，常见的问题是存在电噪声（根据实验室位置不同为 50 或 60 Hz）。当电信号污染不能充分降低时，研究人员会尝试一些方法，如随机化刺激强度和/或使用非整数刺激率来减少其影响。

4. 宽肩概率密度函数

一些生物变量的 pdf 有一个奇特的性质，即距离均值超过三个标准差（SD）的值出现的频率比高斯 pdf 预期的要高得多，也就是观察到的 pdf 有“宽肩”。随着测量设备的改进和数据集的增大，这种现象在自然界中可能是普遍存在的。

这可能会让学过统计学本科课程的读者感到困惑，因为中心极限定理表明，具有有限方差的独立同分布项的归一化和的极限必须是高斯分布。而广义中心极限定理给出了解决这个悖论的方法，它指出独立同分布项的归一化和的唯一可能的非平凡极限是 L´evy - 稳定分布。所有 L´evy - 稳定（或 α - 稳定）分布都有这样的性质：如果 x 和 y 是来自这个家族的随机变量，那么 x + y 形成的随机变量也具有来自这个家族的分布。高斯分布（α = 2）是 α - 稳定分布的特殊情况，此时均值和方差都存在；当 1 < α < 2 时，只有均值存在；当 α ≤ 1 时，均值和方差都不存在。

生物中的随机变量可能由 L´evy 型 pdf 表征，这给实验和建模都带来了挑战。尽管距离均值超过三个标准差的事件发生的可能性增加了，但这些事件仍然非常罕见。因此，必须收集非常大的数据集才能确保检测到这些罕见事件，例如至少需要 10⁵ 个数据点。在收集时间序列时，除非注意低通和高通滤波器的截止频率，否则可能会错过这些罕见事件。而且，采样频率必须足够高，因为采样频率降低时，波动的分布必然会接近高斯分布。

截断是一个特别重要的问题。对于 α 不等于 2 的 L´evy 分布，其尾部延伸到 ±∞，但在现实世界中，物理限制可能会使某些值无法达到，从而截断 L´evy 分布。然而，关于在 L´evy 飞行的背景下如何正确处理截断的研究很少。一种简单的方法是将所有大于预设阈值的 x - x0 设置为零。有研究表明，在面对 L´evy 飞行的控制时，生物会尝试改变截断。

处理 L´evy - 稳定分布的一个实际困难是，只有在三个特殊参数选择下才有 pdf 的封闭形式表达式：高斯分布对应于 α = 2 且 γ = σ²/2 的 L´evy 分布；柯西分布对应于 α = 1 的情况；L´evy 分布对应于 α = 0.5 的情况。相反，L´evy 特征函数 ϕL( f) 对于所有 α 都有封闭形式，因此大多数对 L´evy 飞行性质的研究都使用 ϕL( f)。当 0.5 ≤ α ≤ 2 时，最适合将实验数据拟合到 L´evy 分布的形式为：
[ \phi_{L}(f) = \exp[2\pi j\zeta f - |2\pi\gamma f|^{\alpha} + j\beta g(f,\alpha,\gamma)] ]
其中：
[ g(f,\alpha,\gamma) =
\begin{cases}
2\pi\gamma f \tan\frac{\pi\alpha}{2} (|2\pi f|^{\alpha - 1} - 1) & \text{if } \alpha \neq 1 \
-4\gamma f \log|2\pi\gamma f| & \text{if } \alpha = 1
\end{cases}
]
这里 β 是对称指数，γ 用于缩放分布的范围。

5. 生存函数

生存时间是一个重要的统计量，它在很多领域都有应用，如医学中癌症患者的预期寿命、放射性同位素的寿命、电器或机动车的预期使用寿命以及亚稳态的生存等。

生存时间 tesc 通常被视为随机变量，可以用概率分布来描述。定义累积分布函数 Φ(t) 为状态在时间 t 之前失效的概率：
[Φ(t) = Prob(t_{esc} < t)]
概率密度函数 p(t) 为在小时间间隔内发生失效的概率：
[\varphi(t) = \frac{dΦ(t)}{dt}]
事件持续时间超过时间 t 的概率 p(t) 为：
[p(t) = 1 - Φ(t) = Prob(t_{esc} > t)]
生存函数 s(t) 可以估计 p(t)，通过绘制存活时间超过 t 的个体事件的比例随时间 t 的变化来实验确定。由于 s(t) 近似等于 1 - Φ(t)，因此可以得到 Φ(t) 和 φ(t) 的估计值，进而确定危险函数 h(t)：
[h(t) = \frac{\varphi(t)}{p(t)}]
危险函数给出了单位时间内失效的概率，这个过程被称为 Kaplan - Meier 分析。

通常 p(t) 具有以下形式：
[p(t) = \exp[-(kt)^{\gamma}]]
其中 k 是具有时间倒数单位的常数。当 γ = 1 时，s(t) 近似等于 p(t) = exp(-kt)，φ(t) = kexp(-kt)，h(t) = k，失效风险保持不变；当 γ > 1 时，失效概率随时间增加；当 γ < 1 时，失效概率随时间减小。

例如，在研究指尖平衡杆技能的发展中，生存函数 s(t) 可以通过让受试者进行大量试验（通常 25 - 150 次），并计算杆保持平衡的试验比例随时间的变化来估计。对于新手平衡杆者，在最初 9 天的练习中，技能的提高反映在杆平衡生存曲线逐渐向右移动，即平均杆平衡时间逐渐变长。当时间按 t/t 1/2 重新缩放时，生存曲线会合并为一条曲线，这表明 t 1/2 是杆平衡技能发展的一个相关时间尺度。而且，如果在同一天进行第二次练习，生存曲线没有明显变化，这支持了睡眠依赖的运动控制巩固和细化是技能获取过程的一部分的观点。

6. 相关函数

在时间序列分析中，相关指的是一个变量在给定时间的变化与它过去的行为（自相关函数，用 Cxx(Δ) 表示）或另一个变量的变化（互相关函数，用 Cxy(Δ) 表示）的关联程度。

自相关函数描述了在一个时间测量的变量值可以从之前某个时间测量的值预测的程度。对于确定性系统，其时间演化完全由微分方程描述，不同时间的值是紧密相连的。但如果 x(t) 的值包含随机或随机分量，那么 x(t0) 对 t > t0 时 x(t) 的预测能力就会减弱。因此，测量 Cxx(Δ) 用于估计自我影响变化的时间依赖性，其中滞后 Δ = t - t0。

互相关函数描述了一个变量 x(t) 在时间 t 的值可以从另一个变量 y(t′) 在其他时间 t′ 测量的值预测的程度。因此，测量 Cxy(Δ) 用于估计一个变量的变化与另一个变量的变化的相关程度。需要注意的是，相关并不意味着因果关系。

在分析前，通常会从时间序列中减去均值，所以 Cxx(Δ) 与自协方差函数相同。为了统一，使用“相关”而不是协方差这个术语。用 K 表示自相关或互相关除以 Δ = 0 时的自相关或互相关。

自相关积分 Cxx(Δ) 定义为：
[C_{xx}(\Delta) = \int_{-\infty}^{\infty} x(u)x(u - \Delta)du = \int_{-\infty}^{\infty} x(u)x(u + \Delta)du]
它与卷积积分 z(Δ) 相关，但关键区别在于 Cxx(Δ) 包含 x(u - Δ) 而不是 x(Δ - u)。从图形解释来看，自相关积分没有折叠步骤，所以如果 x(t) 是实函数，Cxx(Δ) 是偶函数，除非 x(t) 本身是偶函数，此时卷积和自相关是等价的。

例如，对于正弦波过程 {xk(t)} = {Asin[2π f0t + φ(k)]}，其中 A、f0 是常数，φ(k) 是从 [0, 2π] 区间上的均匀密度中抽取的随机变量，p(φ) = (2π)⁻¹。计算 Cxx(Δ) 的过程如下：
设
[
\begin{cases}
\alpha = 2\pi f_0t \
\beta = 2\pi f_0(t + \Delta) \
\alpha - \beta = -2\pi f_0\Delta \
\alpha + \beta = 4\pi f_0t + 2\pi f_0\Delta
\end{cases}
]
则
[
\begin{align }
C_{xx}(\Delta) &= A^2 \int_{0}^{2\pi} \sin(\alpha + \varphi)\sin(\beta + \varphi)p(\varphi)d\varphi \
&= \frac{A^2}{2\pi} \left[\frac{2\varphi}{4} \cos(2\pi f_0\Delta) - \frac{\sin(\alpha + \beta + 2\varphi)}{4}\right]_0^{2\pi} \
&= \frac{A^2}{2} \cos(2\pi f_0\Delta)
\end{align }
]

对于实验者来说，自相关随 Δ 的衰减是最有趣的。因此，将 Cxx(Δ) 除以方差 Cxx(0) 得到 Kxx(Δ) 是很有用的。

以下是一个简单的流程图，展示信号平均过程：

graph LR
    A[施加刺激] --> B[记录响应]
    B --> C[多次试验]
    C --> D[奇偶试验平均]
    D --> E[计算信号功率和噪声功率]
    E --> F[提高信噪比]

方法	描述
人工神经网络	用于估计 pdf 的一种方法
核估计	估计 pdf 的方法之一
特征函数法	流行的估计 pdf 方法
最大似然估计器	用于最小化均方误差
信号平均技术	减少噪声影响，增强信号
Kaplan - Meier 分析	用于生存函数分析

随机扰动与相关函数分析

7. 相关函数的性质与应用

相关函数在时间序列分析中具有重要的性质和广泛的应用。自相关函数 (C_{xx}(\Delta)) 是一个偶函数，即 (C_{xx}(\Delta)=C_{xx}(-\Delta))，这意味着变量在时间 (t) 和 (t + \Delta) 的相关性与在时间 (t) 和 (t - \Delta) 的相关性是相同的。这种对称性反映了时间序列在时间上的可逆性（在一定程度上）。

互相关函数 (C_{xy}(\Delta)) 则可以帮助我们发现两个不同变量之间的时间关系。例如，在研究两个不同生理信号（如心电图和脑电图）时，互相关函数可以揭示它们之间是否存在延迟的同步关系。如果 (C_{xy}(\Delta)) 在某个非零的 (\Delta) 值处取得最大值，那么就意味着变量 (y) 的变化相对于变量 (x) 存在一个时间延迟。

在实际应用中，相关函数可以用于以下几个方面：
- 信号检测 ：通过计算自相关函数，可以检测信号中是否存在周期性成分。如果自相关函数在某些特定的 (\Delta) 值处出现明显的峰值，那么就表明信号中存在周期性的波动。
- 系统识别 ：对于一个未知的动态系统，可以通过输入一个已知的信号，然后测量系统的输出，并计算输入和输出之间的互相关函数，从而识别系统的特性。
- 噪声去除 ：如果噪声与信号不相关，那么可以利用相关函数的性质来去除噪声。例如，通过对信号进行多次采样并计算自相关函数，可以减少噪声的影响。

8. 随机扰动下的系统稳定性

随机扰动会对系统的稳定性产生重要影响。在确定性系统中，稳定性通常可以通过分析系统的特征值或李雅普诺夫函数来判断。然而，在存在随机扰动的情况下，系统的稳定性分析变得更加复杂。

考虑一个简单的线性随机系统：
(dx(t) = (Ax(t) + Bw(t))dt)
其中 (x(t)) 是系统的状态向量，(A) 是系统矩阵，(B) 是扰动矩阵，(w(t)) 是随机噪声。

为了分析这个系统的稳定性，我们可以使用随机李雅普诺夫函数。随机李雅普诺夫函数是一个正定函数 (V(x,t))，它满足一定的条件，使得系统的状态在平均意义下是稳定的。具体来说，如果存在一个随机李雅普诺夫函数 (V(x,t))，使得：
(E\left[\frac{dV(x,t)}{dt}\right] \leq 0)
其中 (E[\cdot]) 表示数学期望，那么系统在均方意义下是稳定的。

在实际应用中，随机扰动可能会导致系统出现随机共振现象。随机共振是指在一定的条件下，随机噪声可以增强系统对微弱信号的响应。这种现象在生物系统、物理系统等领域都有广泛的应用。

9. 随机扰动的建模与仿真

为了研究随机扰动对系统的影响，我们需要对随机扰动进行建模和仿真。常见的随机扰动模型包括高斯白噪声、泊松过程等。

高斯白噪声是一种最常见的随机噪声模型，它的特点是在每个时间点上的取值都是独立的，并且服从高斯分布。在 MATLAB 中，可以使用 randn 函数来生成高斯白噪声。

泊松过程是一种用于描述随机事件发生时间的模型。它的特点是事件的发生是独立的，并且在每个时间间隔内发生的事件数服从泊松分布。在 MATLAB 中，可以使用 poissrnd 函数来生成泊松过程。

以下是一个简单的 MATLAB 代码示例，用于模拟一个受高斯白噪声扰动的线性系统：

% 系统参数
A = [-1 0; 0 -2];
B = [1; 1];
T = 10; % 仿真时间
dt = 0.01; % 时间步长
t = 0:dt:T;
N = length(t);

% 初始化状态向量
x = zeros(2, N);
x(:,1) = [1; 1];

% 生成高斯白噪声
w = randn(1, N);

% 仿真系统
for i = 2:N
    dx = (A * x(:,i-1) + B * w(i)) * dt;
    x(:,i) = x(:,i-1) + dx;
end

% 绘制结果
figure;
plot(t, x(1,:), 'b', t, x(2,:), 'r');
xlabel('Time');
ylabel('State');
legend('x_1', 'x_2');

10. 总结与展望

随机扰动在许多领域都有着重要的影响，包括生物学、物理学、工程学等。通过对随机扰动的研究，我们可以更好地理解系统的行为，并开发出更有效的控制和优化策略。

在本文中，我们介绍了随机扰动的相关概念，包括概率密度函数估计、矩的估计、信号平均、宽肩概率密度函数、生存函数和相关函数等。我们还讨论了相关函数的性质和应用，以及随机扰动下的系统稳定性和建模与仿真方法。

未来的研究方向包括：
- 更复杂的随机扰动模型 ：目前的随机扰动模型大多基于简单的假设，如高斯白噪声。未来的研究可以考虑更复杂的随机扰动模型，如非高斯噪声、有色噪声等。
- 随机扰动下的优化问题 ：在存在随机扰动的情况下，如何进行系统的优化是一个重要的问题。未来的研究可以探索随机扰动下的优化算法，如随机梯度下降法、蒙特卡罗树搜索等。
- 随机扰动与深度学习的结合 ：深度学习在许多领域都取得了巨大的成功。未来的研究可以探索随机扰动与深度学习的结合，以提高深度学习模型的鲁棒性和泛化能力。

以下是一个总结表格，列出了本文介绍的主要内容：
| 主题 | 描述 |
| ---- | ---- |
| 概率密度函数估计 | 介绍了估计概率密度函数的方法，包括人工神经网络、核估计和特征函数法 |
| 矩的估计 | 讨论了矩的估计方法，包括无偏估计和有偏估计，以及均方误差的概念 |
| 信号平均 | 介绍了信号平均技术，用于减少噪声影响，增强信号 |
| 宽肩概率密度函数 | 讨论了宽肩概率密度函数的性质和应用，以及 L´evy - 稳定分布 |
| 生存函数 | 介绍了生存函数的概念和应用，以及 Kaplan - Meier 分析方法 |
| 相关函数 | 讨论了自相关函数和互相关函数的性质和应用，以及相关函数在时间序列分析中的作用 |
| 随机扰动下的系统稳定性 | 分析了随机扰动对系统稳定性的影响，以及随机李雅普诺夫函数的应用 |
| 随机扰动的建模与仿真 | 介绍了常见的随机扰动模型，如高斯白噪声、泊松过程等，并给出了 MATLAB 代码示例 |

以下是一个流程图，展示随机扰动研究的一般流程：

graph LR
    A[问题提出] --> B[随机扰动建模]
    B --> C[系统分析]
    C --> D[模型验证与仿真]
    D --> E[结果分析与优化]
    E --> F[应用与推广]