马尔可夫链的模拟答案

【问题标题】：Simulation of Markov chains马尔可夫链的模拟
【发布时间】：2013-11-24 09:28:15
【问题描述】：

我有以下马尔可夫链：

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

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

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

【问题讨论】：

马尔可夫链是基于离散事件的，对吧？除非您有“SimEvents”工具箱，否则这在 Simulink 中可能会变得有点笨重。（或者也许使用这个FEX-alternative，但实际上是用于驱动控制）——无论如何，它可以在 Simulink 中完成，而无需任何进一步的工具箱。这个问题可能很快就会结束，因为它在 SO 上是题外话。

标签： matlab simulink markov-chains

【解决方案1】：

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

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

方法 1：需要计量经济学工具箱
此方法使用dtmc() 和simulate() 函数。

首先，使用概率转移矩阵P 并使用dtmc() 创建Discrete Time Markov Chain (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

【讨论】：

【解决方案2】：

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

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

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

【讨论】：

【解决方案3】：

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

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

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

将其包装在 for 循环中，您可以进行任意数量的模拟。

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

【讨论】：