如何计算两个随机变量的以下概率：R 中的 Pr(Y>X)？

如何解决如何计算两个随机变量的以下概率：R 中的 Pr(Y>X)？

如何在 R 中找到这个概率 P(X

解决方法

cubature 软件包是一组经过试验和测试的用于多变量集成的工具。构建一个包含两个变量的函数作为适当密度函数的参数，并将不等式的条件作为逻辑值包含在内，该逻辑值仅在 X

遗憾的是，有必要为发送到函数的值向量使用单一名称，因此它不像其他方式那样透明：

library(cubature)
prod2beta <- function(x){ (x[1] < x[2]) *    # (X < Y) logical times...
                           dbeta(x[1],1,1) * # X density ...
                           dbeta(x[2],2,3)} # times Y density

hcubature( prod2beta,lower=c(0,0),upper=c(1,1)) # integrate over unit square
#-------------------
 $integral
[1] 0.4

$error
[1] 3.999939e-06

$functionEvaluations
[1] 2168775

$returnCode
[1] 0

这是帮助理解几何情况的线框图：

实现 analytical solution 似乎非常具有挑战性且计算量很大。

如果您对近似解决方案感到满意，请尝试以下任一方法：

方法一：模拟

n <- 1000000
x <- rbeta(n,1)
y <- rbeta(n,3)

sum(x < y)/n
[1] 0.399723

此处的结果取决于您选择的 n 和 RNG，但较高的 n 会产生相当准确的估计。

方法 2：正态近似

方法见Cook (2012)。

x_a <- x_b <- 1
y_a <- 2
y_b <- 3

mu_x <- 1 / (1 + 1)
mu_y <- 2 / (2 + 3)

var_x <- (x_a*x_b) / ( (x_a + x_b)^2 * (x_a + x_b + 1) )
var_y <- (y_a*y_b) / ( (y_a + y_b)^2 * (y_a + y_b + 1) )

pnorm((mu_y - mu_x) / ((var_y + var_x)^0.5))
[1] 0.3879188

这计算量较小，但仍然相当准确。根据库克的说法，这种近似的平均绝对误差为 0.006676，在您的情况下，不高于 0.05069。

我们可以使用 mean 进行蒙特卡罗模拟

> n <- 1e6

> mean(rbeta(n,1) < rbeta(n,3))
[1] 0.400643

正如 Robert Dodier 所指出的，如果所有参数都是整数，那么这个问题在分析上是可以解决的，如果一个分布对于另一个的任何参数都是 Beta(1,1)，那么这个问题就很容易解决。一般来说，如果我们在 [0,1] 上有自变量 X（密度 p，CDF P）和 Y（q 和 Q），那么

P(X<Y) = Integral{ 0<=x<y | p(x)q(y)}
       = Integral{ 0<=y<=1 | Integral{ 0<=x<=y | p(x)} q(y)} 
       = Integral{ 0<=y<=1 | P(y)q(y)}

在手头的情况下 p(x) = 1 所以 P(x) = x

P(X<Y) = Integral{ 0<=y<=1 | y*q(y)}

代替q，

P(X<Y) = Integral( 0<=y<=1 | y*pow(y,a-1)*pow(1-y,b-1)} / B(a,b)
       = Integral( 0<=y<=1 | pow(y,a)*pow(1-y,b-1)}/B(a,b)

但被积函数是 B(a+1,b) 变量的密度，由 B(a+1,b) 缩放，所以

P(X<Y) = B(a+1,b)/B(a,b)

并使用 B 的定义和伽马函数的递归，

P(X<Y) = a/(a+b)

在特定情况下，a=2,b=3,P(X