cubature包是一套久经考验的多变量集成工具。建立一个两个变量的函数作为适当密度函数的参数，并将不等式的条件包括为逻辑值，该逻辑值仅在X < Y的区域中为数值1，并在（X，Y）= {[0，1]，[[0，1]}的联合范围上积分：
遗憾的是，有必要为发送到函数的值的向量使用一个名称，因此它不像其他方式那样透明：

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

这是计算密集度较低的，并且仍然相当准确。根据库克的说法，这种近似的平均绝对误差为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<Y）= 0.4

Beta(1,1)是（0，1）上的均匀分布。我们用U表示。
P(Y > U) = integral P(Y > u) du = integral (1 - F(u)) du其中F是Y的cdf。它是known，这等于E(Y)（Y的期望值），也就是a/(a+b)。