Расширение tecdat|R Language Metropolis Hastings Sampling и байесовская регрессионная модель Пуассона

сбор данных

Оригинальная ссылка:tecdat.cn/?p=23524

Первоисточник:Публичный аккаунт Tuoduan Data Tribe

В этой статье я хочу показать вам, как использовать выборку Metropolis R для выборки из байесовской регрессионной модели Пуассона.

Алгоритм Метрополиса-Гастингса

Алгоритм выборки Метрополиса-Гастингса представляет собой класс методов Монте-Карло с цепями Маркова (MCMC), основная идея которых заключается в создании цепи Маркова.Сделайте его стационарное распределение целевым распределением. Одним из наиболее распространенных применений этого алгоритма является выборка из апостериорных плотностей в байесовской статистике, что и является целью данной статьи.

Алгоритм определяет, как генерировать следующее состояние для данного состояния XtСуществует точка-кандидат Y, которая распространяется из предложения, сгенерированный в , принимается в соответствии с критериями принятия решения, поэтому цепочка переходит в состояние Y в момент времени t+1, т.е. Xt+1=Y, или отвергается, поэтому цепочка остается в состоянии Xt в момент времени t+1, т.е. Xt+ 1=Хт.

Выборка Метрополиса

В алгоритме Метрополиса распределение предложений симметрично, то есть распределение предложенийудовлетворить 

, поэтому процесс сэмплера Метрополиса, генерирующего цепь Маркова, выглядит следующим образом.

  1. Выберите распространение предложения, Прежде чем выбрать его, поймите идеальные характеристики в этой функции.
  2. Сгенерируйте X0 из распределения предложения g.
  3. Повторяйте, пока цепь не сойдется к стационарному распределению.
  • от генерировать Ю.
  • Сгенерируйте U из униформы (0, 1).
  • если , принимает Y и устанавливает Xt+1=Y, иначе устанавливает Xt+1=Xt. Это означает, что точка-кандидат Y принимается с высокой вероятностью.
  • Приращение т.

Байесовский метод

Как я упоминал ранее, мы хотим получить выборку из байесовской модели регрессии Пуассона.

Для оценки параметров в байесовском анализе нам необходимо найти функцию правдоподобия интересующей модели, в данном случае, из модели регрессии Пуассона.

Теперь нам нужно указать априорное распределение для каждого из параметров β0 и β1. Мы будем использовать неинформативные нормальные распределения для этих двух параметров, β0~N(0,100) и β1~N(0,100).

Наконец, мы определяем апостериорное распределение как произведение априорного распределения и распределения правдоподобия.

При использовании пробоотборника Metropolis целевым распределением будет апостериорное распределение.

метод расчета

Здесь вы узнаете, как использовать сэмплер Metropolis R для выборки из апостериорного распределения параметров β0 и β1.

данные

Во-первых, мы генерируем данные из модели регрессии Пуассона, представленной выше.

n <- 1000 #  样本大小
J <- 2 # 参数的数量
X <- runif(n,-2,2) # 生成自变量的值
beta <- runif(J,-2,2) #生成参数的值
y <- rpois(n, lambda = lambda) # 生成因变量的值

Функция правдоподобия

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

LikelihoodFunction <- function(param){
        beta0 <- param[1] 
        beta1 <- param[2] 
        lambda <- exp(beta1*X + beta0)
        # 对数似然函数
        loglikelihoods <- sum(dpois(y, lambda = lambda, log=T)) 
        return(loglikelihoods)
}

предварительное распределение

Далее мы определяем априорные распределения для параметров β0 и β1. Как и в случае с функцией правдоподобия, мы будем использовать журнал априорного распределения.


        beta0prior <- dnorm(beta0, 0, sqrt(100), log=TRUE)
        beta1prior <- dnorm(beta1, 0, sqrt(100), log=TRUE)
        return(beta0prior + beta1prior) #先验分布的对数

Апостериорное распределение

Поскольку мы работаем с логарифмами, мы определяем апостериорное распределение как сумму логарифма функции правдоподобия и логарифма априорного распределения. Помните, что эта функция — наша целевая функция f(.), из которой мы хотим сделать выборку.

Функция предложения

Наконец, мы определяем распределение предложений g(.|Xt). Поскольку мы будем использовать сэмплер Metropolis, распределение предложения должно быть симметричным и зависеть от текущего состояния цепочки, поэтому мы будем использовать нормальное распределение со средним значением, равным значению параметра в текущем состоянии.

Метрополис Сэмплер

Наконец, мы пишем код, который помогает нам запускать сэмплер Metropolis. В этом случае, поскольку мы используем логарифмы, мы должны определить вероятность того, что точка-кандидат Y будет принята как .

        # 创建一个数组来保存链的值
        chain[1, ] <- startvalue # 定义链的起始值
        for (i in 1:iterations){
                # 从提议函数生成Y
                Y <- ProposalFunction(chain[i, ]) 
                # 候选点被接受的概率
                                           PosteriorFunction(chain[i, ]))
                # 接受或拒绝Y的决策标准 
                if (runif(1) < probability) {
                        chain[i+1, ] <- Y
                }else{ 
                        chain[i+1, ] <- chain[i, ]

Из-за сильной автокорреляции цепочки MCMC в краткосрочной перспективе могут быть получены выборки, которые не являются репрезентативными для истинного лежащего в основе апостериорного распределения. Затем, чтобы уменьшить автокорреляцию, мы можем просто разбавить выборку каждым n значениями в цепочке. В этом случае мы будем выбирать значение для нашей последней цепочки каждые 20 итераций алгоритма.

startvalue <- c(0, 0) # 定义链条的起始值
#每20次迭代选择最终链的值
for (i in 1:10000){
        if (i == 1){
                cfinal[i, ] <- chain[i*20,]
        } else {
                cfinal[i, ] <- chain[i*20,]

# 删除链上的前5000个值
burnIn <- 5000

Здесь вы можете увидеть график ACF, который дает нам значение автокорреляции любого ряда со значениями его запаздывания. В этом случае мы показываем график АКФ исходной цепи MCMC и конечной цепи после разбавления образцов для обоих параметров. Из графика можно сделать вывод, что использованная процедура действительно способна значительно уменьшить автокорреляцию.

результат

В этом разделе мы представляем цепочку, созданную пробоотборником Metropolis, и ее распределение по параметрам β0 и β1. Истинное значение параметра представлено красной линией.

Сравнение с glm()

Теперь нам нужно сравнить результаты, полученные с помощью выборки Метрополиса, с функцией glm(), которая используется для подбора обобщенной линейной модели.

В таблице ниже приведены фактические значения параметров и среднее значение оценок, полученных с помощью пробоотборника Metropolis.

##       True value Mean MCMC       glm
## beta0  1.0578047 1.0769213 1.0769789
## beta1  0.8113144 0.8007347 0.8009269

в заключении

Из результатов можно сделать вывод, что оценочные значения параметров β0 и β1 регрессионной модели Пуассона, полученные с помощью пробоотборника Metropolis и функции glm(), очень похожи и близки к фактическим значениям параметров. Кроме того, необходимо признать, что выбор априорного распределения, распределения предложения и начального значения цепочки оказывает большое влияние на результаты, поэтому этот выбор должен быть сделан правильно.

Самые популярные идеи

1.Случай применения многомерной логистической регрессии на языке R

2.Пример реализации панельного регрессионного анализа с плавным переходом (PSTR)

3.Частичная регрессия методом наименьших квадратов (PLSR) и регрессия основных компонентов (PCR) в Matlab

4.R Язычный бассейн Пуассон регрессионный анализ

5.Тест согласия Хосмера-Лемешова в регрессии на языке R

6.Реализация регрессии LASSO, регрессии Ridge и модели эластичной сети на языке r

7.Реализация логистической логистической регрессии на языке R

8.python использует линейную регрессию для прогнозирования цен на акции

9.Как язык R рассчитывает показатели IDI, NRI в анализе выживаемости и регрессии Кокса

Категории