我有以下马尔可夫链:
该链显示了位于小行星带中的宇宙飞船的状态:S1 - 可用,S2 - 已损坏。0.12 - 与小行星碰撞摧毁宇宙飞船的概率。0.88 - 碰撞的概率并不重要。需要计算第三次碰撞后船舶处于可用状态的概率。
分析解决方案显示响应 - 0.681。但是这个问题需要通过任何建模工具(MATLAB Simulink、AnyLogic、Scilab等)的仿真方法来解决。
您知道在 Simulink 或任何其他仿真环境中应该使用哪些组件来仿真这个过程吗?任何示例或链接。
我有以下马尔可夫链:
该链显示了位于小行星带中的宇宙飞船的状态:S1 - 可用,S2 - 已损坏。0.12 - 与小行星碰撞摧毁宇宙飞船的概率。0.88 - 碰撞的概率并不重要。需要计算第三次碰撞后船舶处于可用状态的概率。
分析解决方案显示响应 - 0.681。但是这个问题需要通过任何建模工具(MATLAB Simulink、AnyLogic、Scilab等)的仿真方法来解决。
您知道在 Simulink 或任何其他仿真环境中应该使用哪些组件来仿真这个过程吗?任何示例或链接。
首先,我们知道三步概率转移矩阵包含答案(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
对于更复杂的系统,您需要使用 Stateflow 或 SimEvents,但对于这个简单的示例,您只需要一个带有 Switch 模块的 Unit Delay 模块(输出 = 0 => S1,输出 = 1 => S2),一个随机块和一些比较块来构造确定状态下一个值的逻辑。
据推测,您必须执行(非常)大量的模拟并对结果进行平均以获得具有统计意义的输出。每次运行模拟时,您都需要更改随机发生器的“种子”。这可以通过将种子设置为“现在”(或类似的东西)来完成。
或者,您可以很容易地对模型进行矢量化,这样您只需要执行一次。
如果你想模拟这个,在matlab中是相当容易的:
servicable = 1;
t = 0;
while servicable =1
t = t+1;
servicable = rand()<=0.88
end
现在t
表示船被破坏之前的步数。
将其包装在 for 循环中,您可以根据需要进行尽可能多的模拟。
请注意,这实际上可以为您提供分布,如果您想在 3 次后知道它,只需添加 && t<3
到 while 条件。