MATLAB环境下基于时序蒙特卡罗方法的合成数据生成 基于马尔可夫链蒙特卡罗方法的合成数据生成
最近在帮实验室做时间序列分析的时候,发现用蒙特卡罗方法生成合成数据真是个好用的工具。特别是基于马尔可夫链的这种,特别适合模拟存在状态转移的场景。咱们直接上干货,拿MATLAB举个天气预测的例子——假设天气只有晴天、阴天两种状态,用代码实现状态跳转的过程。
先来定义个靠谱的状态转移矩阵。在MATLAB里可以用二维数组直接表示,每一行代表当前状态,列代表下一状态的概率。比如下面这个矩阵,晴天有80%概率保持晴天,阴天有30%概率转晴:
P = [0.8 0.2; % 晴天→晴天 80%,晴天→阴天 20% 0.3 0.7]; % 阴天→晴天 30%,阴天→阴天 70%接下来是蒙特卡罗模拟的核心部分。这里有个小技巧,用cumsum函数把概率累加起来,再用rand函数生成随机数判断状态跳转。比如这样写循环:
states = zeros(1,1000); % 预分配内存 current_state = randi(2); % 随机初始状态 cumP = cumsum(P,2); % 按行累加概率 for t = 1:1000 r = rand(); if r <= cumP(current_state,1) next_state = 1; else next_state = 2; end states(t) = next_state; current_state = next_state; end这段代码的精髓在于cumsum和rand的配合使用。cumsum把每行的概率累加成[0.8,1.0]和[0.3,1.0]这样的区间,rand生成的随机数掉在哪个区间就对应哪个状态。这样写比用if-else判断更简洁,特别是状态多的时候优势更明显。
生成状态序列后,咱们可以给每个状态绑定观测值。比如晴天对应温度均值为25度,阴天20度,加上随机噪声:
temp = zeros(1,1000); for i = 1:1000 if states(i) == 1 temp(i) = 25 + randn*2; % 晴天温度 else temp(i) = 20 + randn*1.5; % 阴天温度 end end这里有个需要注意的地方:实际应用中要根据领域知识确定噪声的分布类型和参数。比如气温变化可能用正态分布,而股票价格可能用对数正态分布。
验证生成的数据是否合理,可以计算状态转移频率是否接近设定概率。用MATLAB的tabulate函数统计转移次数:
trans_counts = zeros(2); for t = 1:999 i = states(t); j = states(t+1); trans_counts(i,j) = trans_counts(i,j)+1; end disp('实际转移概率:') disp(trans_counts ./ sum(trans_counts,2))运行几次会发现,生成数据的状态转移概率通常在设定值的±5%范围内波动。如果偏差太大,可能需要增加模拟步数或者检查转移矩阵是否满足不可约性条件。
这种生成方法特别适合需要大量训练数据的场景。比如做电力负荷预测时,我见过有人用三层马尔可夫链模拟工作日、周末、节假日的用电模式,每个层级对应不同的状态转移矩阵。只要调整P矩阵的参数,就能快速生成不同季节、不同地区的数据变体。
不过要注意,传统马尔可夫链只能捕捉一阶相关性。如果需要考虑更长期的时间依赖,可以尝试隐马尔可夫模型或者马尔可夫链蒙特卡洛的变体。下次有机会可以聊聊怎么用Metropolis-Hastings算法处理连续状态空间的情况——那又是另一个有意思的故事了。