如何在 R 中找到这个概率 P(X
【问题讨论】：
标签： r probability probability-density probability-distribution
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, 1)
y <- rbeta(n, 2, 3)
sum(x < y)/n
[1] 0.399723
此处的结果取决于您选择的 n 和 RNG，但更高的 n 将产生相当准确的估计。
方法见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
这计算量较小，但仍然相当准确。根据 Cook 的说法，这种近似值的平均绝对误差为 0.006676，在您的情况下，不高于 0.05069。
【讨论】：
我们可以使用mean 进行蒙特卡罗模拟
> n <- 1e6
> mean(rbeta(n, 1, 1) < rbeta(n, 2, 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
【讨论】：