平稳随机过程穿越：从集合平均到单样本概率的降维可视化

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'Arial Unicode MS', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

rng = np.random.default_rng(2023)

# ============================================================
# 生成 100 条独立的窄带平稳随机过程
# ============================================================
fs = 1000.0
T_total = 5.0
t = np.arange(0, T_total, 1/fs)
f0 = 5.0
fc = 0.5
n_samples = 100
a = 0.5  # 阈值

def generate_narrowband(rng):
    n = len(t)
    white = rng.standard_normal(n)
    F = np.fft.rfft(white)
    freqs = np.fft.rfftfreq(n, d=t[1]-t[0])
    F[freqs > fc] = 0
    slow = np.fft.irfft(F, n)
    slow = slow / np.std(slow)
    env = np.abs(slow) + 0.3
    phase = 2*np.pi*f0*t + 0.8*slow
    return env * np.cos(phase)

# 生成集合矩阵：形状为 (100, 5000)
Y = np.array([generate_narrowband(rng) for _ in range(n_samples)])

def find_upcrossings(t, y, a):
    idx = np.where((y[:-1] < a) & (y[1:] >= a))[0]
    return idx

# ============================================================
# 绘图设置：1 行上，2 行下
# ============================================================
fig = plt.figure(figsize=(15, 9))
gs = fig.add_gridspec(2, 2, height_ratios=[1, 1], hspace=0.3, wspace=0.2)

# 选一个特定的观察时刻 t0 和微元 dt
t0 = 2.0
dt = 0.05
idx_t0 = np.argmin(np.abs(t - t0))
idx_dt = np.argmin(np.abs(t - (t0 + dt)))

# ---------- 子图 1：100条曲线的集合视角（宏观） ----------
ax1 = fig.add_subplot(gs[0, :])
# 画前 20 条，颜色调淡
for i in range(20):
    ax1.plot(t, Y[i], color='steelblue', alpha=0.3, lw=0.8)
# 突出显示挑选出的“典型单样本”（第 0 条）
ax1.plot(t, Y[0], color='navy', lw=1.5, label='典型单样本 (Sample 0)')

ax1.axhline(a, color='crimson', ls='--', lw=1.5, label=f'阈值 a={a}')
ax1.axvspan(t0, t0+dt, color='darkorange', alpha=0.3, label=f'微元窗口 dt={dt}s')

ax1.set_xlim(1.5, 3.0)
ax1.set_ylim(-2.5, 2.5)
ax1.set_xlabel(r'时间 $t$ (s)')
ax1.set_ylabel(r'幅值 $y(t)$')
ax1.set_title(r'【宏观视角】100 条集合曲线叠加：在微元 $dt$ 内，我们关注有多少条曲线穿过阈值', fontsize=13)
ax1.legend(loc='upper right')
ax1.grid(alpha=0.3)

# ---------- 子图 2：微元放大图（微观概率） ----------
ax2 = fig.add_subplot(gs[1, 0])
# 截取微元窗口内所有曲线的数据
window_mask = (t >= t0 - 0.02) & (t <= t0 + dt + 0.02)
for i in range(n_samples):
    y_win = Y[i, window_mask]
    t_win = t[window_mask]
    # 找出在这段微元内发生穿越的曲线
    ups = find_upcrossings(t_win, y_win, a)
    if len(ups) > 0:
        ax2.plot(t_win, y_win, color='crimson', alpha=0.6, lw=1.2)
    else:
        ax2.plot(t_win, y_win, color='lightgray', alpha=0.5, lw=0.8)

ax2.axhline(a, color='crimson', ls='--', lw=1.5)
ax2.axvline(t0, color='k', ls=':', lw=1, label='微元起点')
ax2.axvline(t0+dt, color='k', ls=':', lw=1, label='微元终点')
ax2.set_title(r'【微观切片】500 条样本在 $dt$ 内的表现：红线代表发生了穿越 (P = 穿越数/总数)', fontsize=12)
ax2.set_xlabel(r'时间 $t$ (s)')
ax2.set_ylabel(r'幅值 $y(t)$')
ax2.grid(alpha=0.3)

# ---------- 子图 3：单样本概率密度积分（等价三） ----------
ax3 = fig.add_subplot(gs[1, 1])
# 提取所有样本在 y≈a 时的导数 \dot{y}
ydot = np.gradient(Y, t, axis=1)
mask_a = np.abs(Y - a) < 0.05
ydots_at_a = ydot[mask_a]

# 只看正穿越方向的导数 (y_dot > 0)
ydots_pos = ydots_at_a[ydots_at_a > 0]

# 绘制直方图
ax3.hist(ydots_at_a, bins=50, density=True, alpha=0.3, color='gray', label=r'全部斜率分布 $p(\dot{y}|y \approx a)$')
ax3.hist(ydots_pos, bins=30, density=True, alpha=0.7, color='darkorange', label=r'向上穿越部分 $\dot{y} > 0$')

# 标注出选取的典型样本的穿越斜率
# 找到 Sample 0 在窗口内的穿越点
t_win = t[window_mask]
y0_win = Y[0, window_mask]
ups0 = find_upcrossings(t_win, y0_win, a)
if len(ups0) > 0:
    idx_up = ups0[0]
    # 近似计算它的导数
    slope = (y0_win[idx_up+1] - y0_win[idx_up]) / (t_win[idx_up+1] - t_win[idx_up])
    ax3.axvline(slope, color='navy', lw=2, label=fr'Sample 0 的穿越斜率 $\dot{{y}} \approx {slope:.1f}$')

ax3.axvline(0, color='k', lw=1)
ax3.set_xlim(-20, 40)
ax3.set_xlabel(r'导数 $\dot{y}$ (幅值/秒)')
ax3.set_ylabel('概率密度')
ax3.set_title(r'【概率积分】Rice 公式只对 $\dot{y} > 0$ 的正半轴积分', fontsize=12)
ax3.legend(fontsize=9)
ax3.grid(alpha=0.3)

fig.suptitle('平稳随机过程穿越：从集合平均到单样本概率的降维可视化', fontsize=15, y=0.96)
plt.savefig('ensemble_to_single.png', dpi=180, bbox_inches='tight')
plt.show()