条件随机抽样：n个样本必须加到Y

Question

我如何从 X_1 到 X_n（一次性？）对 n 个数字进行采样，使它们的总数是给定的 Y，并且知道 X_i 的无条件分布是（或将是）参数为 mu 和 sigma 的正态分布？

只是绘制 n 个正态分布的数字并将它们缩放到总 Y 似乎不正确。如果 Y 的绝对值非常大，它会倾向于产生 n 个异常值（相对于无条件分布）而不是 - 我不知道 - 1 个异常值。另外：如果 Y 为零，则结果没有意义。

Answer 1

我没试过，但这应该可以用 condMVNorm 包.

让我解释一下 n=2。我们从两个独立的高斯随机变量 X1 和 X2 开始。设 Y=X1+X2。那么很容易得到(X1,Y)的方差-协方差矩阵。然后，使用 condMVNorm，您可以在给定 Y=y 的情况下从 X1 的条件分布中采样。当然，你会得到 X2，X2=Y-X1。

通过设置 Y=X1+...+Xn 并考虑随机向量 (X_1, ...., X_{n-1}, Y)，这可以推广到任何整数 n。

---

编辑

这里是 n=3 的代码。

library(condMVNorm)

# X_i parameters
mu <- 1
sigma <- 2

# (X1, X2, X1+X2+X3) parameters
Mu <- c(mu, mu, 3 * mu)
Sigma <- rbind(
  c(sigma^2, 0,       sigma^2),
  c(0,       sigma^2, sigma^2),
  c(sigma^2, sigma^2, 3*sigma^2)
)

# (X1, X2) given Y=y (where Y = X1+X2+X3)
y <- 5

# number of simulations
nsims <- 10

# simulations of (X1, X2) given Y=y
rcmvnorm(
  nsims, Mu, Sigma, dependent.ind = c(1, 2), given.ind = 3, X.given = y
)

Answer 2

好的，对法线进行采样的有趣方式，但是您必须检查返回的内容

你采样代码

library(condMVNorm)

my_y <- 10
my_sigma <- 5
my_n <- 12

my_mean <- c(rep(my_y/my_n, my_n))
my_bigSigma <- my_sigma^2 * (diag(my_n) - matrix(1, my_n, my_n) / my_n)

my_Result <- rcmvnorm(10000000, my_mean, my_bigSigma, c(1:my_n))

你现在有矩阵

dim(my_Result)

[1] 1000000      12

看起来不错

apply(my_Result, 2, mean)

制作

[1] 0.8373670 0.8266844 0.8293054 0.8366576 0.8358933 0.8316721 0.8415673
[8] 0.8377708 0.8256674 0.8316660 0.8350799 0.8306687

这也可以，但是

apply(my_Result, 2, sd)

会产生

[1] 4.785463 4.786450 4.785305 4.785815 4.781498 4.787159 4.790677 4.788603
[9] 4.788610 4.787577 4.783270 4.786239

和

library(moments)

apply(my_Result, 2, skewness)

会返回类似的东西

[1]  0.0003224842  0.0006144637 -0.0009435844 -0.0016183800  0.0006096481
[6] -0.0016528886 -0.0016519665  0.0020264454 -0.0026967559  0.0009690657
[11]  0.0020904441  0.0016340393

再次看起来不错

apply(my_Result, 2, kurtosis)

会回来

[1] 2.999561 2.999075 3.005832 3.003915 3.006498 3.008577 2.998123 3.000093
[9] 3.008372 3.000561 3.001040 3.000248

接近 3，这对高斯函数有好处

Answer 3

在很多帮助下 (https://stats.stackexchange.com/a/609156/370545)，我认为这可以完成工作：

install.packages("condMVNorm")
library(condMVNorm)

my_y <- 10
my_sigma <- 5
my_n <- 12

my_mean <- c(rep(my_y/my_n, my_n))
my_bigSigma <- my_sigma^2 * (diag(my_n) - matrix(1, my_n, my_n) / my_n)

my_Result <- rcmvnorm(1, my_mean, my_bigSigma, c(1:my_n))[1,]

--

下面的粗略测试，结果sd满足我的要求

install.packages("condMVNorm")
library(condMVNorm)

myout <- c(rep(0, 100000*12))

for (i in 1:100000){
  my_y <- rnorm(1, 10, 5*sqrt(12))
  my_sigma <- 5
  my_n <- 12
  my_mean <- c(rep(my_y/my_n, my_n))
  my_bigSigma <- my_sigma^2 * (diag(my_n) - matrix(1, my_n, my_n) / my_n)
  my_Result <- rcmvnorm(1, my_mean, my_bigSigma, c(1:my_n))[1,]
  myout[((i-1)*12+1):(i*12)] <- my_Result
}

mean(myout)*12
sd(myout)

哪个给

> mean(myout)*12
[1] 9,958756
> sd(myout)
[1] 5,002101

一百万：

> mean(myout)*12
[1] 9,995456
> sd(myout)
[1] 4,999421