r - 如何从真实数据中估计 R 中二元正态分布的参数？

Question

我有一组来自真实数据的 xy 对，我想用二元正态分布建模，由两个正态分布 X 和 Y 组成。我想计算参数，以便我可以重新创建分布而不必使用原始分布源数据，因为它太昂贵(一百万行)。

目前我正在成功地绘制这些数据:

hexbinplot(x~y, data=xyPairs, xbins=16)

我想我需要估计以下参数:

分布平均值 X

分布 X 的标准差

分布平均值 Y

分布 Y 的标准差

Rho，用于创建 Sigma 矩阵

然后双变量法线指定为:



在 R 中有一个包可以做到这一点吗？

我浏览了许多软件包，但其中大多数都可以帮助您使用随机数据模拟双变量，而不是帮助您创建对真实数据进行建模的双变量正态分布。

如果您想了解更多详情，请告诉我。

Answer 1

好的，让我们从几个事实开始:

如果您有一个多元正态分布，边缘分布不依赖于与已被边缘化的变量相关的任何参数。见 here

参数的最大似然估计量 mu和 sigma^2众所周知对应于 sample 类似物。见 here有关如何在单变量情况下获得解析解的示例。

这使我们得出结论，您可以通过以下方式估计这些参数。首先，让我生成一些示例数据:

n <- 10000
set.seed(123) #for reproducible results
dat <- MASS::mvrnorm(n=n, 
               mu=c(5, 10), 
               Sigma= matrix(c(1,0.5,0.5,2), byrow=T, ncol=2)
)

在这里，我选择了 mu1和 mu2分别为 5 和 10。另外， sigma1^2等于 1, rho*sigma1*sigma2等于 0.5，和 sigma2^2等于 2。请注意，由于 rho * sigma1 * sigma2 = 0.5 ，我们有 rho = 0.5/sqrt(1*2) = 0.35
使用已知(分析)最大似然估计

现在，让我们估计参数 mu1和 mu2首先从数据。在这里，我使用每个单独变量的样本均值，因为事实 1 确保我不需要担心依赖性。也就是说，我可以忽略它们是双变量正态的，因为边缘分布具有相同的参数，而且我碰巧知道在单变量情况下这些参数的 MLE 是样本均值。

> colMeans(dat)
[1]  5.006143 9.993642

我们看到这非常接近我们之前在生成数据时指定的真实值。

现在，让我们估计 x1 的方差和 x2 :

> apply(dat, 2, var)
[1] 0.9956085 2.0008649

此外，这非常接近真实值。到目前为止，这种方法似乎效果很好。 :)

现在，剩下的就是 rho :注意方差协方差矩阵的非对角线上的条目是 rho*sigma1*sigma2 = rho * 1 * sqrt(2) ，我定义为 0.5。因此， rho = 0.35 .

现在，让我们看一下样本相关性。样本相关已经标准化了协方差，所以我们不需要手动除以 sqrt(2)得到相关系数。

 > cor(dat)
          [,1]      [,2]
[1,] 1.0000000 0.3481344
[2,] 0.3481344 1.0000000

这再次非常接近先前指定的 true 参数。请注意，有人可能会争辩说后者在小样本中存在偏差，我们可以进行更正。有关讨论，请参阅维基百科文章。如果你想这样做，你只需将最后一项乘以 n/(n-1) .样本大小如 n=10000 ，它通常不会产生很大的不同。

现在，我在这里做了什么？我知道这些量的分析最大似然估计量是什么样子的，我只是用它们来估计这些参数。如果您不知道解决方案在分析上的样子，您会怎么做？原则上，您知道似然函数。你有数据。您可以将似然函数写为参数的函数，然后只需使用许多可用优化器中的一种来查找使样本似然最大化的参数值。这将是直接的机器学习方法。见 here .

所以，让我们试试吧。

在数值上最大化似然

上述过程利用了我们能够通过分析获得最大似然估计量的事实。也就是说，我们通过取似然函数的导数，将其设为零，并求解未知量，找到了这些量的封闭形式解。但是，我们也可以使用计算机以数值方式查找值，如果您找不到易于处理的解析解，这可能会派上用场。让我们试试看。

首先，由于我们要最大化一个函数，让我们使用内置函数 optim为了那个原因。 optim要求我提供一个带有初始起始值的参数向量，以及一个将参数向量作为参数的函数。该函数应该返回一个要最大化或最小化的值。

该函数将是样本似然。给定一个大小为 n 的 iid 样本，样本似然是所有 n 的乘积个体可能性(即概率密度函数)。大产品的数值优化是可能的，但人们通常采用对数将产品转化为总和。要获得似然性，只需仔细观察二元正态分布的个体 pdf，您就会看到样本似然性可以写为

-n*(log(sig1) + log(sig2) + 0.5*log(1-rho^2)) - 
0.5/(1-rho^2)*( sum((x1-mu1)^2)/sig1^2 + 
                sum((x2-mu2)^2)/sig2^2 - 
                2*rho*sum((x1-mu1)*(x2-mu2))/(sig1*sig2) )

此函数将在其参数上最大化。自 optim要求我提供一个参数向量，我为此使用了一个包装器并将最大化问题设置如下:

wrap <- function(parms, dat){
  mymu1  = parms[1]
  mymu2  = parms[2]
  mysig1 = parms[3]
  mysig2 = parms[4]
  myrho  = parms[5]
  myx1 <- dat[,1]
  myx2 <- dat[,2]
  n = length(myx1)

  f <- function(x1=myx1, x2=myx2, mu1=mymu1, mu2=mymu2, sig1=mysig1, sig2=mysig2, rho=myrho){
    -n*(log(sig1) + log(sig2) + 0.5*log(1-rho^2)) - 0.5/(1-rho^2)*(
      sum((x1-mu1)^2)/sig1^2 + sum((x2-mu2)^2)/sig2^2 - 2*rho*sum((x1-mu1)*(x2-mu2))/(sig1*sig2)
      )
  }
  f(x1=myx1, x2=myx2, mu1=mymu1, mu2=mymu2, sig1=mysig1, sig2=mysig2, rho=myrho)

  }

我的电话 optim然后看起来如下:

eps <- eps <- .Machine$double.eps  # get a small value for bounding the paramter space to avoid things such as log(0).

numML <- optim(rep(0.5,5), wrap, dat=dat, 
         method="L-BFGS-B", 
         lower = c(-Inf, -Inf, eps, eps, -1+eps), 
         upper = c(Inf, Inf, 100, 100, 1-eps), 
         control = list(fnscale=-1))

在这里， rep(0.5,5)提供起始值， wrap在函数之上， lower和 upper是参数的界限， fnscale参数确保我们最大化函数。结果，我得到:

numML$par 
[1]  5.0061398 9.9936433 0.9977539 1.4144453 0.3481296

请注意，这些元素对应于 mu1 , mu2 , sig1 , sig2和 rho .如果你方 sig1和 sig2 ，您会看到我们重新创建了我最初提供的差异。所以，它似乎有效。 :)