Синтаксис для "распространяется как" в R - PullRequest
Новая
       15

Синтаксис для "распространяется как" в R

1 голос
/ 10 апреля 2020

При указании модели в JAGS / BUGS символ «распространяется как» ~ очень полезен. Как это сделать в R при использовании методов MCM C, которые требуют от меня указать вероятность?

Допустим, я хочу оценить три параметра, которые являются многомерными нормально распределенными. В JAGS я бы сделал это, указав pars[1:n] ~ dmnorm(mu[1:3], sigma[1:3, 1:3]). Если все указано правильно, JAGS включит go для оценки этих параметров при заданном распределении.

В R есть аналогичные функции, например, функция dmvnorm() из пакета mvtnorm. Однако я не уверен, как их использовать. Я должен предоставить данные, чтобы получить плотность вероятности, тогда как в JAGS мне нужно только предоставить параметры распределения, такие как mu и sigma. Что такое R-эквивалент синтаксиса ~ в JAGS?

1 Ответ

1 голос
/ 10 апреля 2020

Вот некоторые случайные данные: 

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 это бета (1,1), которая эквивалентна равномерному распределению.

Так что давайте (первоначально) забудем байесовскую часть. Вы спросили, как рассчитать вероятность. Вероятность для параметра тета заданных независимых и одинаково распределенных данных y = (y_1, ..., y_N) составляет

L (theta | y) = product (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
>

Мы можем заставить его работать для всего диапазона значений p, используя sapply

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

Так что теперь, например, я мог вычислить вероятность значений p в диапазоне от 0,01 до 0,99 с шагом 0,01 на 

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

И я мог бы построить их 

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

enter image description here

Вы можете увидеть из Сюжет, что вероятность максимальна при 0,3, поэтому это MLE p на основе этих данных.

Возвращаясь к проблеме Байеса, вот реализация Метрополис-Гастингс (без комментариев, извините): 


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)
}

Теперь мы возьмем образец размером 10000 из этого, с * 104 8 * 

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 образца (плюс несколько достоверных интервалов)

enter image description here

И увидите, что Метрополис-Гастингс имеет сработало - интервалы велики, потому что размер выборки невелик.

...