乘概率分布函数

Multiply Probability Distribution Functions

我很难构建一个有效的程序,通过对概率密度函数进行加法和乘法来预测完成两个过程步骤所需的时间分布。

让"a"表示完成过程"A"需要多长时间的概率分布函数。零天 = 10%,一天 = 40%,两天 = 50%。令 "b" 表示完成过程 "B" 所需时间的概率分布函数。零日 = 10%,一日 = 20%,等等

进程 "B" 在进程 "A" 完成之前无法启动,因此 "B" 依赖于 "A"。

a <- c(.1, .4, .5)
b <- c(.1,.2,.3,.3,.1)

如何计算完成"A"和"B"的时间概率密度函数?

这是我期望的输出或以下示例:

totallength <- 0 # initialize
totallength[1:(length(a) + length(b))] <- 0 # initialize
totallength[1] <- a[1]*b[1]
totallength[2] <- a[1]*b[2] + a[2]*b[1]
totallength[3] <- a[1]*b[3] + a[2]*b[2] + a[3]*b[1]
totallength[4] <- a[1]*b[4] + a[2]*b[3] + a[3]*b[2]
totallength[5] <- a[1]*b[5] + a[2]*b[4] + a[3]*b[3]
totallength[6] <- a[2]*b[5] + a[3]*b[4]
totallength[7] <- a[3]*b[5]

print(totallength)
[1] [1] 0.01 0.06 0.16 0.25 0.28 0.19 0.05
sum(totallength)
[1] 1

我在 Visual Basic 中有一种使用三个 for 循环的方法(一个用于每个步骤,一个用于输出)但我希望我不必在 R 中循环。

因为这似乎是一个非常标准的流程问题,我的问题的第二部分是是否存在任何库来模拟操作流程,所以我不会从头开始创建它。

执行此类操作的有效方法是使用卷积:

convolve(a, rev(b), type="open")
# [1] 0.01 0.06 0.16 0.25 0.28 0.19 0.05

这是高效的,因为它比单独计算每个值的输入更少,还因为它以高效的方式实现(使用快速傅里叶变换或 FFT)。

您可以使用您发布的公式确认这些值中的每一个都是正确的:

(expected <- c(a[1]*b[1], a[1]*b[2] + a[2]*b[1], a[1]*b[3] + a[2]*b[2] + a[3]*b[1], a[1]*b[4] + a[2]*b[3] + a[3]*b[2], a[1]*b[5] + a[2]*b[4] + a[3]*b[3], a[2]*b[5] + a[3]*b[4], a[3]*b[5]))
# [1] 0.01 0.06 0.16 0.25 0.28 0.19 0.05

我不熟悉完全符合您的示例描述的专用包。但让我为这个问题提供一个更强大的解决方案。 您正在寻找一种方法来估计可能由 n 个步骤过程组合的过程的分布,在您的案例 2 中可能不像您的示例那样容易计算。 我将使用的方法是模拟,从底层分布中淹没 10k 个观测值,然后计算模拟结果的密度函数。 使用您的示例,我们可以执行以下操作:

x <- runif(10000)
y <- runif(10000)

library(data.table)
z <- as.data.table(cbind(x,y))
z[x>=0 & x<0.1, a_days:=0]
z[x>=0.1 & x<0.5, a_days:=1]
z[x>=0.5 & x<=1, a_days:=2]
z[y>=0 & y <0.1, b_days:=0]
z[x>=0.1 & x<0.3, b_days:=1]
z[x>=0.3 & x<0.5, b_days:=2]
z[x>=0.5 & x<0.8, b_days:=3]
z[x>=0.8 & x<=1, b_days:=4]
z[,total_days:=a_days+b_days]
hist(z[,total_days])

如果你的第二个进程被指数分布淹没,密度和方法也可以工作,这将产生一个很好的代理。在这种情况下,您将使用 rexp 函数直接计算 b_days。

查看包裹:distr。选择术语 "multiply" 是不幸的,因为所描述的情况不是对概率的贡献是独立的(其中概率的乘法将是自然使用的术语)。它更像是某种顺序加法,而这正是 distr 包提供的,它解释了“+”在用作两个离散分布的符号操作时的含义。

 A <- DiscreteDistribution ( setNames(0:2, c('Zero', 'one', 'two') ), a)
 B <- DiscreteDistribution(setNames(0:2, c(  "Zero2" ,"one2", "two2", 
                                               "three2", "four2") ),  b )
?'operators-methods'  # where operations on 2 DiscreteDistribution are convolution
plot(A+B)

经过一番搜索后,我发现实际的数值可以在这里找到:

 A.then.B <- A + B
> environment(A.the.nB@d)$dx
[1] 0.01 0.06 0.16 0.25 0.28 0.19 0.05

似乎应该有一种显示概率的方法,而且我不是这个迷人软件包的普通用户,所以很可能有一个。一定要阅读小插图和代码演示……我还没有读完。进一步四处闲逛使我确信正确的地方是在配套包装中:distrDoc,其中的插图有 100 多页长。而且它也不应该需要任何努力来找到它,因为该建议在加载包时打印的消息中......除了在我的辩护中有几页消息,所以它更诱人跳转到编码和使用帮助页面。