如何用蒙特卡洛方法生成同维度含NaN的非平稳时序伪随机值

我有一个一维非平稳（含周期性+趋势）时间序列Ts，其中包含NaN值。需要基于该序列的概率分布，生成10000个与Ts同维度（保留原序列中NaN与非NaN的位置）的伪随机值。

部分数据 #

import numpy as np
NaN = np.nan
Ts = np.array([384.540,378.233,376.858,378.497,NaN,NaN,NaN,NaN,NaN,NaN,NaN,390.409,386.174,382.2768,382.082,383.721,NaN,NaN,NaN,NaN,NaN,NaN,NaN,391.841,389.513,382.835,381.387,384.404,NaN,NaN,NaN,NaN,NaN,NaN,NaN,393.871,391.176,385.041,385.270,385.570,NaN,NaN,NaN,NaN,NaN,NaN,NaN,398.377,395.187,390.173,387.628,388.129,NaN,NaN,NaN,NaN,NaN,NaN,NaN,395.886,390.830,389.398,391.617,NaN,NaN,NaN,NaN,NaN,NaN,399.943,390.400,391.019,393.635,NaN,NaN,NaN,NaN,NaN,NaN,403.128,399.594,394.948,394.561,395.420,NaN,NaN,NaN,NaN,NaN,NaN,NaN,405.345,403.449,398.429,395.195,397.791,NaN])

已尝试的代码 #

rr=(Ts-np.nanmean(Ts))/np.nanstd(Ts); # 归一化
mu, sigma = np.nanmean(rr), np.nanstd(rr) # 均值和标准差
q=np.random.uniform(mu, sigma, rr.shape[0]); # 生成均匀分布随机值

原序列是非平稳的（含趋势和周期性），直接用全局统计量生成的随机值无法匹配原序列的分布特征。基于蒙特卡洛思路，正确的做法是依托原序列非NaN值的经验分布抽样，同时严格保留原序列的NaN位置，具体实现步骤如下：

步骤1：提取有效数据并构建经验分布 #

从原序列中提取所有非NaN的有效值，基于这些值构建经验分布作为抽样依据：

import numpy as np
from scipy.stats import rv_histogram

# 提取非NaN有效值
valid_vals = Ts[~np.isnan(Ts)]

# 基于有效值构建经验分布
hist, bins = np.histogram(valid_vals, bins='auto', density=True)
empirical_dist = rv_histogram((hist, bins))

步骤2：批量生成10000个同维度样本 #

对每个样本，保留原序列的NaN位置，在非NaN位置填充从经验分布中抽取的随机值：

# 获取原序列的NaN位置掩码
nan_mask = np.isnan(Ts)
n_samples = 10000
# 初始化样本数组，形状为(10000, 原序列长度)
samples = np.full((n_samples, len(Ts)), np.nan)

# 批量填充抽样值
for i in range(n_samples):
    # 生成与非NaN数量匹配的抽样值
    sampled_vals = empirical_dist.rvs(size=len(valid_vals))
    # 将抽样值填充到对应非NaN位置
    samples[i, ~nan_mask] = sampled_vals

进阶优化（若需保留时序相关性） #

如果需要保留原序列的时序依赖特征（比如趋势、周期性），仅匹配边缘分布不够，可改用：

• 时序自举法：对原序列的非NaN连续时序块进行随机抽样重组
• 时序模型拟合：先用ARIMA、Prophet等模型拟合原序列的趋势和周期，对残差进行蒙特卡洛抽样后，再还原趋势和周期得到完整样本

结果验证 #

可以通过对比原序列与生成样本的分布验证效果：

import matplotlib.pyplot as plt

# 绘制原序列有效值分布
plt.hist(valid_vals, bins='auto', density=True, alpha=0.5, label='原序列')
# 绘制单个生成样本的非NaN值分布
plt.hist(samples[0, ~nan_mask], bins='auto', density=True, alpha=0.5, label='生成样本')
plt.legend()
plt.show()

内容的提问来源于stack exchange，提问作者Shubho