如何使用R将纸上定义的统计模型“翻译”到计算机上？

编程问答 2022-08-10

问题描述

我最初将此问题发布在stats.stackexchange.com上， 但由于专注于编程，因此已关闭。希望我 可以在这里获得任何帮助。

为了简化起见，我不会在此处提供许多理论上的细节，但我的最终目标是使用R实现隐马尔可夫模型。

尽管我对理论模型的构建很满意，但是当我尝试实现它时，我意识到我不了解计算统计的基本知识。我的问题就是朝这个方向发展。

让

$X$

和

$Y$

是随机变量，例如

enter image description here

和

$Y|X \sim \mathcal{N}(x,\sigma^2)$

，以及

$p \in [0,1]$

和

$\sigma^2 \in [0,+\infty)$

。如果

$\pi(\cdot)$

表示分布，则如何计算

enter image description here

使用 R？

我的意思是，这些分布（一离散和一连续）乘法的确切含义是什么？如何使用R执行此操作？答案显然是

$x$

的函数，但是它在我的代码中如何表示？

如果

$Y|X$

也是离散的，是否有任何变化？例如，

enter image description here

和

$q \in [0,1]$

。它将如何影响已实现的代码？

我知道我的问题不是很具体，但是我对如何开始非常迷失。我对这个问题的目标是了解如何将纸上写的内容“翻译”到计算机上。

解决方法

翻译

等式描述了在观察到X以及参数Y=y和p的值的情况下如何计算sigma的概率分布。最终，您想要实现一个函数p_X_given_Y，该函数采用值Y并返回X的概率分布。一个好的开始是实现表达式的RHS中使用的两个功能。像

p_X <- function (x,p=0.5) { switch(as.character(x),"0"=p,"1"=1-p,0) }

p_Y_given_X <- function (y,x,sigma=1) { dnorm(y,sd=sigma) }

请注意，此处任意选择了p和sigma。然后可以使用这些函数来定义p_X_given_Y函数：

p_X_given_Y <- function (y) {
  # numerators: for each x \in X
  ps <- sapply(c("0"=0,"1"=1),function (x) { p_X(x) * p_Y_given_X(y,x) })

  # divide out denominator
  ps / sum(ps)
}

可以这样使用：

> p_X_given_Y(y=0)
#         0         1 
# 0.6224593 0.3775407

> p_X_given_Y(y=0.5)
#   0   1 
# 0.5 0.5 

> p_X_given_Y(y=2)
#         0         1 
# 0.1824255 0.8175745

这些数字应具有直观意义（给定p=0.5）：Y=0更有可能来自X=0，Y=0.5也有可能X=0或X=1等。这只是实现它的一种方法，其思想是返回“ X的分布”，在这种情况下，它只是一个命名的数字矢量，其中名称（“ 0” ，“ 1”）对应于X的支持，而值对应于概率质量。

一些替代实现可能是：

  • 一个p_X_given_Y(x,y)，它也取x的值并返回相应的概率质量
  • 一个p_X_given_Y(y)，该函数返回另一个接受x自变量并返回相应概率质量的函数（即概率质量函数）
bayesianhidden-markov-modelsmachine-learningrr