机器算法验证 - 是否有双变量ββ分布我可以适合我的数据吗？ - 吾爱随笔录

是否有双变量ββ分布我可以适合我的数据吗？

机器算法验证 r 密度函数贝塔分布系词双变量

2022-03-15 00:11:25

我正在分析二维数据。在 fitdistrplus 和 logspline 包的帮助下分析每个维度后，它们都符合 Beta 分布。是否可以像双变量 Beta 分布一样分析二维数据？（注意：我使用的是 R。）
数据集样本：
数据点是随着时间对特定产品进行的 2 种不同化学反应测试的输出。所以在时间 1 PB=2.394和DBA=134.417，在时间 2 PB=2.594和DBA=111.382等等。

structure(list(PB = c(2.394, 2.594, 2.78, 2.499, 2.478, 2.744, 
2.563, 2.553, 2.631, 2.434, 2.604, 2.685, 2.439, 2.548, 2.778, 
2.604, 2.638, 2.585, 2.808, 2.784, 2.489, 2.797, 2.619, 2.687, 
2.624, 2.341, 2.712, 2.493, 2.616, 2.562), DBA = c(134.417, 111.382, 
125.303, 163.445, 89.428, 141.881, 140.559, 141.408, 122.498, 
128.099, 115.88, 111.83, 170.685, 89.956, 128.948, 131.064, 170.114, 
101.843, 132.092, 173.86, 91.976, 130.882, 132.016, 157.143, 
122.052, 100.08, 140.079, 144.167, 141.072, 143.787)), .Names = c("PB", 
"DBA"), row.names = c(NA, 30L), class = "data.frame")

Beta 分布的缩放样本数据集：

    structure(list(PB = c(0.589027911453321, 0.781520692974013, 0.960538979788258, 
0.690086621751685, 0.669874879692012, 0.925890279114534, 0.751684311838306, 
0.742059672762271, 0.817131857555342, 0.627526467757459, 0.791145332050048, 
0.869104908565929, 0.632338787295477, 0.737247353224254, 0.958614051973051, 
0.791145332050048, 0.823869104908566, 0.772858517805582, 0.987487969201155, 
0.964388835418672, 0.68046198267565, 0.976900866217517, 0.8055822906641, 
0.871029836381136, 0.810394610202118, 0.538017324350337, 0.895091434071223, 
0.684311838306064, 0.80269489894129, 0.750721847930703), NOH = c(0.371624288211084, 
0.241754524440435, 0.320240175903479, 0.535282178496927, 0.117979365168856, 
0.413705812707899, 0.406252466595253, 0.411039070868805, 0.304425776625134, 
0.336003833793764, 0.267113942605852, 0.244280317979365, 0.576100806224277, 
0.120956193268309, 0.340790438067317, 0.352720302193156, 0.572881547048543, 
0.18797429103005, 0.358516096295879, 0.594001240345042, 0.13234481592152, 
0.351694198567965, 0.358087613463382, 0.499751930991712, 0.301911258950217, 
0.17803461690252, 0.403546259232114, 0.426594125274849, 0.409144725714608, 
0.424451711112364)), .Names = c("PB", "NOH"), row.names = c(NA, 
30L), class = "data.frame")

1个回答

定义二元贝塔分布的方法有很多种，即平方上的二元分布 $[0,1]\times [0,1]$ 使用 beta 边际。一种方法是从使用 gamma 变量的 beta 分布的通常随机表示开始，让 $X\sim\mathcal{Gamma}(\alpha,\theta),~~ Y\sim\mathcal{Gamma}(\beta,\theta)$ （独立），那么

\frac{X}{X + Y} \sim B e t a (α, β)

$\frac{X}{X+Y} \sim\mathcal{Beta}(\alpha,\beta)$ 如果我们可以用三个 gamma 变量进行类似的表示，并在两个比率中使用其中一个，我们将得到一个带有 beta 边际的相关双变量分布。但请注意，这样的结构不能给出负相关！一篇讨论这个和许多其他结构的论文（arXiv）给出了以下表示

X = \frac{G_{1}}{G_{0} + G_{1}}, Y = \frac{G_{2}}{G_{0} + G_{2}}

$X=\frac{G_1}{G_0+G_1},\quad Y=\frac{G_2}{G_0+G_2}$ 带形状参数

α_{i}

$\alpha_i$ 和常用的尺度参数。密度函数为：

f (x, y) = \frac{1}{B (α_{0}, α_{1}, α_{2})} \frac{x^{α_{1} - 1} (1 - x)^{α_{0} + α_{2} - 1} y^{α_{2} - 1} (1 - y)^{α_{0} + α_{1} - 1}}{(1 - x y)^{α_{0} + α_{1} + α_{2}}}

$f(x,y)=\frac{1}{B(\alpha_0,\alpha_1,\alpha_2)}\frac{x^{\alpha_1-1}(1-x)^{\alpha_0+\alpha_2-1}y^{\alpha_2-1}(1-y)^{\alpha_0+\alpha_1-1}}{(1-xy)^{\alpha_0+\alpha_1+\alpha_2}}$ 在哪里

B

$B$ 是一个广义的 beta 函数

B (α_{1}, α_{2}, α_{3}, \dots) = \frac{\prod_{j} Γ (α_{j})}{Γ (\prod_{j} α_{j})}

$B(\alpha_1,\alpha_2,\alpha_3, \dots)=\frac{\prod_j \Gamma(\alpha_j)}{\Gamma(\prod_j \alpha_j)}$ 我们可以使用最大似然估计对您的数据进行测试，该数据具有较低的正相关性。下面的一些R代码：

gbeta <- function(a, log=FALSE) { # Generalized beta function
    p <- length(a); if(p<2) stop("a must have length at least 2")
    if (min(a)<=0)stop("negative parameters is not allowed")
    lgbeta <- sum(lgamma(a)) - lgamma(sum(a))
    if(log)lgbeta else exp(lgbeta)
    }

dbibeta3 <- function(x, y, a0, a1, a2, log=FALSE) {
    logdens <- -gbeta(c(a0, a1, a2), log=TRUE)  +  (a1-1)*log(x) +
        (a0 + a2-1)*log(1-x) + (a2-1)*log(y) + (a0 + a1-1)*log(1-y) -
        (a0 + a1 + a2)*log(1-x*y)
    if(log) logdens else exp(logdens)
}

rbibeta3 <- function(n, a0, a1, a2) {
    if(min(c(a0, a1, a2))<= 0) stop("all parameters must be positive")
    G0 <- rgamma(n, a0, 1)
    G1 <- rgamma(n, a1, 1)
    G2 <- rgamma(n, a2, 1)
    cbind(G1/(G1 + G0), G2/(G2 + G0))
}

betadata <- < data from your question >

library(bbmle) # on CRAN

minusloglik <- function(a0, a1, a2) { # we use log parameters
    x <- betadata[, 1]; y <- betadata[, 2]
    -sum(dbibeta3(x, y, exp(a0), exp(a1), exp(a2), log=TRUE))
    }

mod <- bbmle::mle2(minusloglik, start=list(a0=log(1), a1=log(3), a2=log(7)))
> mod

Call:
bbmle::mle2(minuslogl = minusloglik, start = list(a0 = log(1), 
    a1 = log(3), a2 = log(7)))

Coefficients:
       a0        a1        a2 
0.7116677 2.0963634 0.2249110 

Log-likelihood: 27.22 
> exp(coef(mod))
      a0       a1       a2 
2.037386 8.136527 1.252211

让我们用这些参数模拟一些数据，并将它们与您的数据一起绘制：

simdata <- rbibeta3(30, 2.037, 8.137, 1.252)

模拟数据显示出比观察到的更强的相关性，所以也许这个模型不是一个好的模型？我们可以尝试一些更灵活的双变量 beta 分布！

其它你可能感兴趣的问题

上一篇Q 学习与价值迭代下一篇如何确定我的深度学习模型需要什么类型的层？