Выборка из настраиваемой функции правдоподобия

У меня есть следующая функция правдоподобия, которую я использовал в довольно сложной модели (на практике в логарифмическом масштабе):

library(plyr)
dcustom=function(x,sd,L,R){
    R. = (log(R) - log(x))/sd
    L. = (log(L) - log(x))/sd
    ll = pnorm(R.) - pnorm(L.)
    return(ll) 
}

df=data.frame(Range=seq(100,500),sd=rep(0.1,401),L=200,U=400)
df=mutate(df, Likelihood = dcustom(Range, sd,L,U))

with(df,plot(Range,Likelihood,type='l'))
abline(v=200)
abline(v=400)

В этой функции sd предопределено, а L и R являются «наблюдениями» (очень похожими на конечные точки равномерного распределения), поэтому указаны все 3 из них. Вышеупомянутая функция обеспечивает большую вероятность (1), если оценка модели x (производный параметр) находится между диапазоном LR, плавное уменьшение правдоподобия (от 0 до 1) около границ (из которых резкость зависит от sd) , и 0, если снаружи слишком много.

Эта функция очень хорошо работает для получения оценок x, но теперь я хотел бы сделать обратное: нарисовать случайный x из приведенной выше функции. Если бы я сделал это много раз, я бы сгенерировал гистограмму, которая повторяет форму кривой, построенной выше.

Конечная цель - сделать это на C ++, но я думаю, что для меня было бы проще, если бы я сначала понял, как это сделать на R.

В Интернете есть полезная информация, которая помогает мне начать работу (http://matlabtricks.com/post-44/generate-random-numbers-with-a-given-distribution, https://stats.stackexchange.com/questions/88697/sample-from-a-custom-continuous-distribution-in-r), но я Я все еще не совсем уверен, как это делать и как это кодировать.

Я предполагаю (совсем не уверен!) Шаги следующие:

  1. преобразовать функцию правдоподобия в распределение вероятностей
  2. вычислить кумулятивную функцию распределения
  3. выборка с обратным преобразованием

Это правильно, и если да, то как мне это закодировать? Спасибо.


r log-likelihood
person Wave    schedule 26.02.2018    source источник

Ответы (1)


arrow_upward
2
arrow_downward

Одна из идей может заключаться в использовании алгоритма Metropolis Hasting для получения выборки из распределения с учетом всех других параметров и вашей вероятности.

# metropolis hasting algorithm
 set.seed(2018)
 n_sample <- 100000
 posterior_sample <- rep(NA, n_sample)
 x <- 300 # starting value: I chose 300 based on your likelihood plot
 for (i in 1:n_sample){
     lik <- dcustom(x = x, sd = 0.1, L = 200, R =400)
     # propose a value for x (you can adjust the stepsize with the sd)
     x.proposed <-  x + rnorm(1, 0, sd = 20)
     lik.proposed <- dcustom(x = x.proposed, sd = 0.1, L = 200, R = 400)
     r <- lik.proposed/lik # this is the acceptance ratio
         # accept new value with probablity of ratio 
         if (runif(1) < r) {
         x <- x.proposed
     posterior_sample[i] <- x
     }
    }

 # plotting the density  
 approximate_distr <- na.omit(posterior_sample) 
 d <- density(approximate_distr)
 plot(d, main = "Sample from distribution")
 abline(v=200)
 abline(v=400)

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

# If you now want to sample just a few values (for example, 5) you could use 
 sample(approximate_distr,5)
#[1] 281.7310 371.2317 378.0504 342.5199 412.3302
person Daniel    schedule 26.02.2018
comment
Здорово!! Спасибо. Есть оптимальный sd (размер шага)? - person Wave; 27.02.2018
comment
Какой размер шага выбрать, всегда зависит от конкретного случая. Проще говоря, если размер шага установлен слишком большим, слишком много значений будет отклонено в процедуре, если размер шага слишком мал, у вас есть опасность, что вы не охватите весь дистрибутив (например, застряли на локальном максимуме). Я бы посоветовал дополнительно прочитать о диагностике mcmc и поиграть с параметрами настройки как для этого случая, так и для других более известных распределений (например, распределения Гаусса), чтобы увидеть, что будут делать различные значения параметров настройки. :) - person Daniel; 27.02.2018