Как я могу “перевести” статистическую модель, определенную на бумаге, на компьютер с использованием R?

#r #machine-learning #bayesian #hidden-markov-models

#r #машинное обучение #байесовский #скрытые марковские модели

Вопрос:

Я изначально разместил этот вопрос на stats.stackexchange.com , но он был закрыт из-за сосредоточенности на программировании. Надеюсь, я смогу получить любую помощь здесь.

Я не буду приводить здесь много теоретических деталей, чтобы упростить ее, но моя конечная цель — реализовать скрытую марковскую модель с использованием R .

Хотя я в порядке с построением теоретической модели, когда я попытался ее реализовать, я понял, что не знаю элементарных вещей о вычислительной статистике. Мой вопрос касается этого направления.

Пусть $X $ и $Y$ будут случайными величинами, такими, что введите описание изображения здесь и $ Y | X  sim  mathcal{N}(x,  sigma ^ 2) $, с $ p  in [0, 1] $ и $  sigma ^ 2  in [0,    infty) $. Если $  pi( cdot) $ обозначает распределение, как я могу вычислить

введите описание изображения здесь

используя R ?

Я имею в виду, каково точное значение умножения этих распределений (одно дискретное и одно непрерывное)? Как я могу это сделать с помощью R ? Ответ, очевидно, зависит от $x $, но как он представлен в моем коде?

Есть ли какие-либо изменения, если $Y|X$ также дискретно? Например, введите описание изображения здесь, с $ q  in [0, 1]$ помощью. Как это повлияет на реализованный код?

Я знаю, что мои вопросы не очень конкретны, но я очень запутался, с чего начать. Моя цель в этом вопросе — понять, как я могу «перевести» то, что я написал на бумаге, на компьютер.

Ответ №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, x, 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 аргумент и возвращает соответствующую массу вероятности (т. Е. Функцию массы вероятности)