从多元学生 t 分布中提取

机器算法验证 r 随机生成

2022-03-20 14:39:08

所以我想弄清楚从多元学生t分布中采样是否有一个很好的分解，就像从多元正态分布中采样一样： http ://en.wikipedia.org/wiki/Multivariate_normal_distribution#Drawing_values_from_the_distribution

有很多R包可以做到这一点，但对我来说，它们都不是很直观的实现。鉴于我具有指定数量的自由度、已知的均值和协方差矩阵，我试图从多元 t 分布中进行采样。

我需要能够从多元 t 分布中采样，因为我已经从我的后验分布中分析整合了参数。即积分后得到。由于我对有一个封闭的分析形式，我希望能够直接从多元学生 t 分布中采样。

p (ψ | y) = \int_{R^{p}} \int_{R^{+}} p (β, σ, ψ | y) d σ^{2} d β

$p(\psi|y)=\int_{R^p}\int_{R^+}p(\beta,\sigma,\psi|y)d\sigma^2d\beta$

ψ \sim Multivariate Student-t (μ, Σ, ν)

$\psi\sim\text{Multivariate Student-t}(\mu,\Sigma,\nu)$

ψ

$\psi$

2个回答

关于多元 t的维基百科文章解释说，当具有分布并且独立地具有分布时，则 $Y$ $N(0, \Sigma)$ $U$ $\chi_\nu^2$

X = μ + \frac{Y}{\sqrt{\frac{U}{ν}}} = μ + Y \sqrt{\frac{ν}{U}}

$X = \mu + \frac{Y}{\sqrt{\frac{U}{\nu}}} = \mu + Y\sqrt{\frac{\nu}{U}}$

具有分布。 $t_\nu(\mu, \Sigma)$

这概括了学生 t 的单变量情况，它是标准正态和（缩放的）变量的比率：我们在这个公式中认识到偏差（预期为零）与标准误差的经典比率，即由学生 t 检验的原假设预测。 $\chi$

( -vectors) 的实现可以从 (也是 -vectors) 和 (numbers) 的独立实现的相同代数组合中获得，因此代价是为每个的系数：这相当有效。 $X$ $p$ $Y$ $p$ $U$ $p+1$ $p$ $X$

此R代码说明了的情况。 $p=2$

library(MASS) # mvrnorm()
mu <- c(7, 5); sigma <- matrix(c(1, 1/2, 1/2, 1), 2); nu <- 4
n <- 4e3 # Number of draws
y <- t(t(mvrnorm(n, rep(0, length(mu)), sigma) * sqrt(nu / rchisq(n, nu))) + mu)

您可以检查colMeans(y)近似值mu，var(y)近似值sigma乘以（Stéphane Laurent 提供的一个观察结果，他的说法是合理的，因为Student t 分布的方差是），你可以看到分布。 $\nu/(\nu-2)$ $\nu/(\nu-2)$ plot(y)

使用四种不同的包和绘图进行模拟：

# variance-covariance matrix :
Sigma <- matrix(rbind(c(2,1),
                      c(1,3)),
                nrow=2)
# mean 
Mu <- c(0,0)
# degrees of freedom
df <- 30
# number of simulations
nsims <- 1000
# sampling function 
sampling <- function(pkg){
  switch(pkg, 
    "mnormt" = mnormt::rmt(nsims, mean=Mu, S=Sigma, df=df),
    "mvtnorm" = mvtnorm::rmvt(nsims, sigma=Sigma, df=df) +
      matrix(rep(Mu,nsims), ncol=2, byrow=TRUE), 
    "LaplacesDemon" = LaplacesDemon::rmvt(nsims, mu=Mu, S=Sigma, df=df),
    "mixAK" = mixAK::rMVT(nsims, df=df, mu=Mu, Sigma=Sigma)$x
  )
}

# sampling
pkgs <- c("mnormt", "mvtnorm", "LaplacesDemon", "mixAK")
sims <- vapply(pkgs, 
               FUN=sampling, FUN.VALUE=matrix(0, nrow=nsims, ncol=2))

# dataframe for ggplot2
dat <- setNames(cbind(rep(pkgs,each=nsims), 
                      data.frame(do.call(rbind, plyr::alply(sims, 3)))),
                c("pkg", "x", "y"))

# plot
library(ggplot2)
ggplot(dat, aes(x=x, y=y)) +
  geom_point() +
    facet_grid(.~pkg)

笔记：

mvtnorm::rmvt允许非中心参数向量。

基准

> library(microbenchmark)
> microbenchmark(
+   mnormt = sampling("mnormt"),
+   mvtnorm = sampling("mvtnorm"),
+   LaplacesDemon = sampling("LaplacesDemon"),
+   mixAK = sampling("mixAK"),
+   times=100
+ )
Unit: microseconds
          expr     min      lq     mean   median       uq      max neval cld
        mnormt 327.989 346.509 402.5122 360.5655 383.9930 3537.376   100 a  
       mvtnorm 518.982 554.235 615.9864 571.4160 604.8840 3814.493   100  b 
 LaplacesDemon 375.737 408.536 461.7511 422.5930 437.5425 3646.705   100 a  
         mixAK 875.530 911.677 938.2057 927.0720 948.2685 1237.434   100   c

更新

我刚刚发现了这个mvnfast包裹，它是赢家。

df <- 5
p <- 6
Mu <- rep(1, p)
Sigma <- toeplitz(p:1)
sampling <- function(pkg, nsims=1000){
  switch(pkg, 
         "mnormt" = mnormt::rmt(nsims, mean=Mu, S=Sigma, df=df), 
         "mvtnorm" = sweep(mvtnorm::rmvt(nsims, sigma=Sigma, df=df), 2, Mu, "+"), 
         "mvnfast" = mvnfast::rmvt(nsims, mu=Mu, sigma=Sigma, df=df)
  )
}

> microbenchmark(
+   mnormt = sampling("mnormt"),
+   mvtnorm = sampling("mvtnorm"),
+   mvnfast = sampling("mvnfast"),
+   times=1000
+ )
Unit: microseconds
         expr      min       lq      mean    median       uq       max neval
       mnormt  757.723  858.573 1156.3124  908.5530 1377.555  5285.759  1000
      mvtnorm 1240.112 1408.346 1956.2496 1527.7160 2253.086 98021.043  1000
      mvnfast  493.099  528.799  599.8516  554.4585  673.828  3693.114  1000

其它你可能感兴趣的问题

上一篇使用核函数是否使数据线性可分？如果是这样，为什么要使用软边距 SVM？下一篇二元变量的预测