机器算法验证 - 泊松/负二项式后有或无替换抽样的解析求解 - 吾爱随笔录

泊松/负二项式后有或无替换抽样的解析求解

机器算法验证数理统计泊松分布联合分配近似

2022-04-05 11:22:28

精简版

我正在尝试分析解决/近似由独立泊松绘制和进一步抽样（有或没有替换）产生的复合可能性（我真的不在乎哪一个）。我想使用 MCMC (Stan) 的可能性，所以我只需要一个常数项的解决方案。最终，我想模拟一个初始抽签来自 neg 的过程。二项式分布，但我想我可以通过泊松情况的解决方案到达那里。

解决方案很可能不可行（我对数学的理解不足以判断这是一个简单的问题还是非常困难的问题）。因此，我也对近似值、负面结果或直觉为什么问题可能难以处理（例如，与已知的难题相比）感兴趣。可以帮助我前进的有用论文/定理/技巧的链接是很好的答案，即使它们与手头问题的联系没有完全解决。

正式声明

更正式地说，首先是独立绘制的，然后我中随机个项目得到。即我从一个瓮中画出个彩色球，其中颜色的球的数量是从中提取的。在这里，已知且固定，我们以为条件。从技术上讲，采样是在没有更换的情况下完成的，但假设有更换的采样应该没什么大不了的。 $Y = (y_1, ..., y_N), y_n \sim Pois(\lambda_n)$ $k$ $Y$ $Z = (z_1,...,z_N)$ $k$ $n$ $Pois(\lambda_n)$ $k$ $\sum_n y_n \geq k$

我尝试了两种方法来解决没有替换的采样问题（因为这似乎是更容易的情况，因为某些术语被取消了），但都被两种方法都卡住了。无放回抽样的可能性是：

P (Z = (z_{1}, . . ., z_{N}) | Λ = (λ_{1}, . . ., λ_{N})) = \frac{\sum_{Y; \forall n : y_{n} \geq z_{n}} (\frac{\prod_{n = 1}^{N} (\binom{y_{n}}{z_{n}})}{(\binom{\sum_{n = 1}^{N} y_{n}}{k})} \prod_{n = 1}^{N} P o i s s o n (y_{n} | λ_{n}))}{P (\sum_{n} y_{n} \geq k | Λ)}

$P(Z = (z_1, ..., z_N) | \Lambda = (\lambda_1, ..., \lambda_N)) = \frac{ \sum_{Y;\forall n: y_n \geq z_n} \left( \frac{\prod_{n=1}^N{y_n \choose z_n}}{ {\sum_{n=1}^N y_n} \choose k} \prod_{n=1}^N Poisson(y_n |\lambda_n) \right) }{ P(\sum_n y_n \geq k|\Lambda) }$

编辑： “尝试的解决方案”部分已删除，因为答案中的解决方案没有建立在它们之上（并且更好）

1个回答

无需更换的情况

如果你有独立的泊松分布变量 $n$

Y_{i} \sim Pois (λ_{i})

$Y_i \sim \text{Pois$(\lambda_i)$}$

和条件

\sum_{j = i}^{n} Y_{j} = K

$\sum_{j=i}^n Y_j = K$

然后

{Y_{i}} \sim Multinom (K, (\frac{λ_{i}}{\sum_{j = 1}^{n} λ_{j}}))

$\lbrace Y_i \rbrace \sim \text{Multinom} \left(K,\left(\frac{\lambda_i}{\sum_{j=1}^n \lambda_j} \right)\right)$

次彩色球填充您的骨灰盒，例如首先绘制总的值（这是条件的泊松分布截止），然后根据多项分布 $n$ $Y_i$ $K$ $K \geq k$ $K$

个球填充瓮相当于为每个球独立绘制分类分布的颜色。然后，您可以将已添加到瓮中（当该样本被绘制而没有替换时）并且其分布只是另一个多项式分布向量： $K$ $k$ $\lbrace Z_i \rbrace$

{Z_{i}} \sim Multinom (k, (\frac{λ_{i}}{\sum_{j = 1}^{n} λ_{j}}))

$\lbrace Z_i \rbrace \sim \text{Multinom} \left(k,\left(\frac{\lambda_i}{\sum_{j=1}^n \lambda_j} \right)\right)$

模拟

##### simulating sample process for 3 variables #######


# settings
set.seed(1)
k = 10
lambda = c(4, 4, 4)
trials = 1000000

# observed counts table
Ocounts = array(0, dim = c(k+1, k+1, k+1))

for (i in c(1:trials)) {
  # draw poisson with limit sum(n) >= k
  repeat {
    Y = rpois(3,lambda)
    if (sum(Y) >= k) {break}
  }
  # setup urn
  urn <- c(rep(1, Y[1]), rep(2, Y[2]), rep(3, Y[3]))
  # draw from urn
  s <- sample(urn, k, replace=0)
  Z = c(sum(s==1),sum(s==2),sum(s==3))
  Ocounts[Z[1]+1, Z[2]+1, Z[3]+1] = Ocounts[Z[1]+1, Z[2]+1, Z[3]+1] + 1
}



# comparison
observed = rep(0, 0.5*k*(k+1))
expected = rep(0, 0.5*k*(k+1))   
t = 1

for (z1 in c(0:k)) {
  for (z2 in c(0:(k-z1))) {  
    z3 = k-z1-z2 
    observed[t] = Ocounts[z1+1, z2+1, z3+1]
    expected[t] = trials*dmultinom(c(z1, z2, z3), prob=lambda)
    t = t+1
  }
}

plot(observed,expected)
x2 <- sum((observed-expected)^2/expected)
pvalue <- 1-pchisq(x2, 66-1)

结果

> # results from chi-sq test
> x2
[1] 75.49286
> pvalue
[1] 0.1754805

负二项式

对于负二项分布（在某些条件下）导致狄利克雷多项分布的情况，这些论点的作用相同。

下面是一个示例模拟

##### simulating sample process for 3 variables #######

# dirichlet multinomial for vectors of size 3
ddirmultinom =  function(x1,x2,x3,p1,p2,p3) {
  (factorial(x1+x2+x3)*gamma(p1+p2+p3)/gamma(x1+x2+x3+p1+p2+p3))/
  (factorial(x1)*gamma(p1)/gamma(x1+p1))/
  (factorial(x2)*gamma(p2)/gamma(x2+p2))/
  (factorial(x3)*gamma(p3)/gamma(x3+p3))
}

# settings
set.seed(1)
k = 10
theta = 1
lambda = c(4,4,4)
trials = 1000000

# calculating negative binomials pars
means = lambda
vars = lambda*(1+theta)

ps = (vars-means)/(vars)
rs = means^2/(vars-means)


# observed counts table
Ocounts = array(0, dim = c(k+1,k+1,k+1))

for (i in c(1:trials)) {
  # draw poisson with limit sum(n) >= k
  repeat {
    Y = rnbinom(3,rs,ps)
    if (sum(Y) >= k) {break}
  }
  # setup urn
  urn <- c(rep(1,Y[1]),rep(2,Y[2]),rep(3,Y[3]))
  # draw from urn
  s <- sample(urn,k,replace=0)
  Z = c(sum(s==1),sum(s==2),sum(s==3))
  Ocounts[Z[1]+1,Z[2]+1,Z[3]+1] = Ocounts[Z[1]+1,Z[2]+1,Z[3]+1] + 1
}



# comparison
observed = rep(0,0.5*k*(k+1))
expected = rep(0,0.5*k*(k+1))   
t = 1

for (z1 in c(0:k)) {
  for (z2 in c(0:(k-z1))) {  
    z3 = k-z1-z2 
    observed[t]=Ocounts[z1+1,z2+1,z3+1]
    expected[t]=trials*ddirmultinom(z1,z2,z3,lambda[1]/theta,lambda[2]/theta,lambda[3]/theta)
    t = t+1
  }
}

plot(observed,expected)
x2 <- sum((observed-expected)^2/expected)
pvalue <- 1-pchisq(x2,66-1)

# results from chi-sq test
x2
pvalue

其它你可能感兴趣的问题