自制蒙特卡洛模拟在安全风险分析中的应用（第二部分）

之前我写过如何将Doug Hubbard的《如何衡量网络安全中的一切》书中的简单定量分析实现为JavaScript代码。当我编写蒙特卡洛模拟代码时，我使用的是基于预期发生率的百分比概率，这在发生率较高时（分析期内发生两次或更多次）显得不够灵活。

去年晚些时候，我与Doug简短交流后，他建议我使用泊松分布，或者缩短分析时间范围。后者看起来像是一种临时解决方案，尤其是在同一投资组合中，有些风险是每五年发生一次，而有些则预期每天发生一次。

由于Excel中没有现成的泊松分布逆采样功能，我转而使用R编程语言。

我从R中的qpois函数开始，该函数能够对分位数进行采样。以下调用根据预期发生率（lambda）在每次采样点提供发生率。runif()函数生成一个0到1之间的随机数，用于选择从泊松分布的哪个分位数进行采样：

1
2


occurrence <- 36
qpois(runif(1,0:1),lambda=occurrence)

R提供了一个名为replicate()的强大函数，允许将多个采样结果存入列表（重要的是每次采样都会重新生成随机数）。以下代码片段对泊松分布进行了1000次采样，并创建了一个包含100个样本的列表：

1
2
3


occurrence <- 36
sample <- 1000
occurrencelist <- replicate(sample, qpois(runif(1,0:1),lambda=occurrence), simplify = TRUE)

这是一个相对简单的蒙特卡洛模拟的一半。现在我们需要考虑风险发生后的后果。之前我借鉴了Doug使用对数正态分布来估计损失的方法，该分布有一些研究支持。然而，对数正态分布具有长尾，可能产生一些较大的数值，特别是在结果范围不确定性较高时，可能会导致意外的结果。

修正的PERT分布允许我们在考虑不确定估计的同时，保留对数正态分布的偏斜特性，但对长尾值的处理更为谨慎。在之前的实现中，我实际上已经对长尾值进行了截断。mc2d包中的qpert()函数提供了一个分位数采样，可以像泊松采样一样进行随机化。

1
2
3
4
5


minharm <- 10000
likelyharm <- 50000
maxharm <- 250000
confidence <- 4
qpert(runif(1,0:1),min=minharm, mode=likelyharm, max=maxharm, shape=confidence)

将泊松分布的发生率采样输出与修正PERT分布的一组采样结果结合，可以按以下方式进行：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21


sample <- 10000
occurrence <- 36
minharm <- 10000
likelyharm <- 50000
maxharm <- 250000
confidence <- 4

occurrencelist <- replicate(sample, qpois(runif(1,0:1),lambda=occurrence), simplify = TRUE)

outcomelist <- vector("numeric", length(occurrencelist))
x <- 1
for (i in occurrencelist)
{
  loopvector <- unlist(replicate(i, qpert(runif(1,0:1),min=minharm, mode=likelyharm, max=maxharm, shape=confidence),simplify=FALSE))

  average <- sum(loopvector)/length(loopvector)
  outcomelist[x] <- average
  x = x+1
}

print(mean(na.omit(outcomelist)))

这段代码运行了10000次模拟采样，并打印了结果的均值（平均值）。虽然代码简单粗糙，但它是有效的。希望这能帮助其他走类似道路的人。Netflix发布了一个名为RiskQuant的成熟开源项目，如果你在寻找工具而不是自己编写代码，我推荐使用它。

一如既往，如果你发现任何明显的错误，或者只是想告诉我我的代码有多糟糕，你可以给我留言（顺便说一句，我知道我的代码有多糟糕）。