如何将贝叶斯计算视为 R 中的干净矩阵

How to see Bayesian computations as clean matrices in R

鉴于理解 R 如何帮助进行贝叶斯计算所需的清晰度,在下文中,我将在这方面提出 R 编码问题。

一些必要的细节:

假设,我有一个名为 mu 的对象,定义为:

mu <- rnorm( 1e4 , 178 , 20 )         ## A vector of hypothesized values

对象 mu 将作为下一个名为 y.given.mu 的对象的 mean 参数:

y.given.mu <- rnorm( 1e4 , mu , 1 )   ## A vector of normal densities conditional on `mu`

问题

我想知道我怎么能:

A) 看清y.given.mu?

的矩阵结构

B)将对象mu乘以y.given.mu,可以清楚地看到这两个对象的乘积矩阵结构(即联合分布)

C)B) 中整合出 mu 这样我就得到了 p(y)?

正如我们讨论的那样,我们将您上一个问题的所有后续问题 移到另一个线程中。


合理足够的网格

delta.mu <- 0.5    # this affects numerical integration precision
mu <- seq(178 - 3 * 20, 178 + 3 * 20, by = delta.mu)
delta.y <- 1    # this does not affect precision, by only plotting
y <- seq(min(mu) - 3, max(mu) + 3, by = delta.y)

# the range above is chosen using 3-sigma rule of normal distribution.
# normal distribution has near 0 density outside (3 * sd) range of its mean

条件密度p(y | mu)

cond <- outer(y, mu, dnorm)
dimnames(cond) <- list(y = y, mu = mu)

# each column is a conditional density, conditioned on some `mu`
# you can view them by for example `plot(y, cond[, 1]), type = "l")
# you can view all of them by `matplot(y, cond, type = "l", lty = 2)`

关节密度p(y, mu)

# marginal of `mu`
p.mu <- dnorm(mu, 178, 20)
# multiply `p.mu` to `cond` column by column (i.e., column scaling)
joint <- cond * rep(p.mu, each = length(y))

边际密度p(y)

# numerical integration by Simpson / Trapezoidal Rule
p.y <- rowSums(joint * delta.mu)

现在让我们绘制并检查

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