matlab - 马尔可夫链的模拟

Question

我有以下马尔可夫链：

该链显示了位于小行星带中的宇宙飞船的状态：S1 - 可用，S2 - 已损坏。0.12 - 与小行星碰撞摧毁宇宙飞船的概率。0.88 - 碰撞的概率并不重要。需要计算第三次碰撞后船舶处于可用状态的概率。

分析解决方案显示响应 - 0.681。但是这个问题需要通过任何建模工具（MATLAB Simulink、AnyLogic、Scilab等）的仿真方法来解决。

您知道在 Simulink 或任何其他仿真环境中应该使用哪些组件来仿真这个过程吗？任何示例或链接。

Answer 1

首先，我们知道三步概率转移矩阵包含答案（0.6815）。

% MATLAB R2019a
P = [0.88 0.12;
    0 1];
P3 = P*P*P
P(1,1)                 % 0.6815

方法 1：需要 Econometrics Toolbox
此方法使用dtmc()和simulate()函数。

首先，创建具有概率转移矩阵的离散时间马尔可夫链 (DTMC)P，并使用dtmc()。

mc = dtmc(P);               % Create the DTMC
numSteps = 3;               % Number of collisions

您可以使用 轻松获得一个示例路径simulate()。注意你如何指定初始条件。

% One Sample Path
rng(8675309)                                     % for reproducibility
X = simulate(mc,numSteps,'X0',[1 0])

% Multiple Sample Paths
numSamplePaths = 3;
X = simulate(mc,numSteps,'X0',[numSamplePaths 0])  % returns a 4 x 3 matrix

第一行是 DTMC 的起始状态（初始条件）的 X0 行。第二行是 1 次转换 (X1) 后的状态。因此，第四行是 3 次转换（碰撞）后的状态。

% 50000 Sample Paths
rng(8675309)                                     % for reproducibility
k = 50000;
X = simulate(mc,numSteps,'X0',[k 0]);            % returns a 4 x 50000 matrix

prob_survive_3collisions = sum(X(end,:)==1)/k    % 0.6800

我们可以在 3 次碰撞中幸存的平均概率上引导 95% 的置信区间以获得0.6814 ± 0.00069221，或者更确切地说，[0.6807 0.6821]，其中包含结果。

numTrials = 40;
ProbSurvive_3collisions = zeros(numTrials,1);
for trial = 1:numTrials
    Xtrial = simulate(mc,numSteps,'X0',[k 0]);
    ProbSurvive_3collisions(trial) = sum(Xtrial(end,:)==1)/k;
end

% Mean +/- Halfwidth
alpha = 0.05;
mean_prob_survive_3collisions = mean(ProbSurvive_3collisions)
hw = tinv(1-(0.5*alpha), numTrials-1)*(std(ProbSurvive_3collisions)/sqrt(numTrials))
ci95 = [mean_prob_survive_3collisions-hw mean_prob_survive_3collisions+hw]   

---

maxNumCollisions = 10;
numSamplePaths = 50000;
ProbSurvive = zeros(maxNumCollisions,1);
for numCollisions = 1:maxNumCollisions
    Xc = simulate(mc,numCollisions,'X0',[numSamplePaths 0]);
    ProbSurvive(numCollisions) = sum(Xc(end,:)==1)/numSamplePaths;
end

Answer 2

对于更复杂的系统，您需要使用 Stateflow 或 SimEvents，但对于这个简单的示例，您只需要一个带有 Switch 模块的 Unit Delay 模块（输出 = 0 => S1，输出 = 1 => S2），一个随机块和一些比较块来构造确定状态下一个值的逻辑。

据推测，您必须执行（非常）大量的模拟并对结果进行平均以获得具有统计意义的输出。每次运行模拟时，您都需要更改随机发生器的“种子”。这可以通过将种子设置为“现在”（或类似的东西）来完成。

或者，您可以很容易地对模型进行矢量化，这样您只需要执行一次。

Answer 3

如果你想模拟这个，在matlab中是相当容易的：

servicable = 1;
t = 0;
while servicable =1
   t = t+1;
  servicable = rand()<=0.88
end

现在t表示船被破坏之前的步数。

将其包装在 for 循环中，您可以根据需要进行尽可能多的模拟。

---

请注意，这实际上可以为您提供分布，如果您想在 3 次后知道它，只需添加 && t<3到 while 条件。