R 中 "is distributed as" 的语法

Question

在JAGS/BUGS中指定模型时，"is distributed as"符号~非常有用。当使用需要我指定可能性的 MCMC 方法时，如何在 R 中执行此操作？

比方说，我想估计三个呈多元正态分布的参数。在 JAGS 中，我将通过指定 pars[1:n] ~ dmnorm(mu[1:3], sigma[1:3, 1:3]) 来执行此操作。如果一切都正确指定，JAGS 将继续在给定分布下估计这些参数。

在 R 中，有类似的函数，例如来自 mvtnorm 包的 dmvnorm() 函数。但是，我不确定如何使用这些。我必须提供数据以获得概率密度，而在 JAGS 中，我只需要提供分布参数，如 mu 和 sigma。与 JAGS 中的 ~ 语法等效的 R 是什么？

Answer 1

这是一些随机数据：

set.seed(123)
y = rbinom(10, 1, 0.2)
y

> y
 [1] 0 0 0 1 1 0 0 1 0 0

因此我们知道生成此数据的 p 值为 0.2。让我们看看如何尝试恢复该信息（假设我们不知道）。在 JAGS 中，我会编写以下模型：

model{
  for(i in 1:10){
    y[i] ~ dbern(p)
  }
  p ~ dunif(0, 1)
}

所以我说过数据是使用具有参数 p 的伯努利分布（或从中采样）生成的，并且 p 的先验是 Beta(1,1)，它相当于均匀分布。

所以让我们（最初）忘记贝叶斯部分。您已经询问了如何计算可能性。给定独立同分布数据 y = (y_1, ..., y_N) 的参数 theta 的似然度是

L(theta | y) = 乘积(f(y_i | theta), i = 1,...,N)

在我们的示例中，pdf f(y_i | theta) 是 p^y_i * (1 - p)^(1 - y_i)。我知道这只是简化为 p 如果 y_i 为 1，或 (1 - p) 如果 y_i 为零，但假设我们不知道这一点并且我们只是使用二项式概率函数参数 n = 1 和 p 来计算这个，那么你可以得到这样的可能性：

Like = function(p){
  prod(dbinom(y, 1, p))
}

这是一个非常简单的函数，仅适用于 p 的单个值，但它确实有效，例如

> Like(0.1)
[1] 0.0004782969
> Like(0.2)
[1] 0.001677722
> Like(0.3)
[1] 0.002223566
>

我们可以使用 sapply

使其适用于 p 的整个值范围

Like = function(p){
  sapply(p, function(p.i)prod(dbinom(y, 1, p.i)))
}

所以现在，例如，我可以通过

以 0.01 的步长计算 p 值范围从 0.01 到 0.99 的可能性

p = seq(0.01, 0.99, by = 0.01)
l = Like(p)

我可以绘制它们

plot(p, l, type = "l")

从图中可以看出，似然在 0.3 处最大化，因此这是基于此数据的 p 的 MLE。

回到贝叶斯问题，这是 Metropolis-Hastings 的实现（抱歉，未注释）：


MH = function(N = 1000, p0 = runif(1)){
  log.like = function(p){
    sum(dbinom(y, size = 1, p, log = TRUE))
  }

  ll0 = log.like(p0)
  r = c(p0, rep(0, N))

  for(i in 1:N){
    p1 = runif(1)
    ll1 = log.like(p1)

    if(ll1 > ll0 || log(runif(1)) < ll1 - ll0){
      p0 = p1
      ll0 = ll1
    }

    r[i + 1] = p0
  }

  return(r)
}

现在我们从中抽取一个大小为 10,000 的样本，

set.seed(123)
p = MH(10000)  
plot(density(p))    
abline(v = c(mean(p), mean(p) + c(-1,1)*qnorm(0.975)*sd(p)))

并绘制样本的 KDE（加上一些可信区间）

并且看到 Metropolis-Hastings 起作用了——间隔很宽，因为样本量很小。

R 中 "is distributed as" 的语法

Syntax for "is distributed as" in R

r

distribution

jags