1

我有以下马尔可夫链:

在此处输入图像描述

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

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

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

4

3 回答 3

2

首先,我们知道三步概率转移矩阵包含答案(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
于 2020-11-03T04:51:02.773 回答
1

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

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

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

于 2013-11-25T06:06:42.067 回答
0

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

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

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

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


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

于 2013-12-24T11:06:04.740 回答