Я хочу сделать выборку из апостериорной, где LambdaA и LambdaB - экспоненциальные скорости A и B. Кроме того, y - это наблюдения случайных величин.
Апостериор дается
и по числовым причинам я беру журнал этой функции.
Данные:
n<-100
y<- c(rexp(n))
Логарифм апостериорного:
dmix<-function(LambdaA,LambdaB,w){
ifelse( LambdaA<=0|LambdaB<=0|w<0|w>1 ,0,log(w*LambdaA*LambdaB*exp(-2*(LambdaA+LambdaB))*prod(w*LambdaA*exp(-
LambdaA*y) + (1-w)*LambdaB*exp(-LambdaB*y)) ))}
U-значения
U.lambdaB <- runif(1)
U.lambdaA<- runif(1)
U.w<- runif(1)
Считать шаги
REJLambdaB <- 1
REJw <- 1
REJLambdaA<-1
Начальные точки
LambdaB <- LambdaA<- w<- numeric(n)
LambdaA[1]<-0.5
LambdaB[1] <- 0.5
w[1] <- 0.5
Алгоритм случайного блуждания MH, обновляющий каждый компонент за раз:
for (t in 2:n){
LambdaBprop<- rnorm(1,LB[t-1],0.5)
wprop<- rnorm(1,w[t-1],0.5)
LambdaAprop<- rnorm(1,LB[t-1],0.5)
logalpha1 = dmix(LambdaAprop,LambdaB[t-1],w[t-1])-dmix(LambdaA[t-1],LambdaB[t-
1],w[t-1])
logalpha2 = dmix(LambdaA[t-1],LambdaBprop,w[t-1])-dmix(LA[t-1],LB[t-1],w[t-
1])
if (!is.null(log(U.lambdaB) > logalpha2))
{LambdaB[t] <- LambdaBprop} ## accepted
else{LambdaB[t] <- LambdaB[t-1] ##rejected
REJLambdaB<-REJLambdaB+1}
if (!is.null(log(U.lambdaA) > logalpha1))
{LambdaA[t]<-LambdaAprop}
else {LambdaA[t]<-LambdaA[t-1]
REJLambdaA<-REJLambdaA+1}
if (w[t]<0|w[t]>1)
{w[t]<-w[t-1]}
else {w[t]<-wprop
REJw<-REJw+1}
}
В конечном счете, у меня проблемы с апостериором, так как я продолжаю получать либо бесконечность, либо 0 при оценке логарифма. Обратите внимание, что я хочу сравнить log($\alpha(x'|x))$ с log(U). Любая помощь, чтобы заставить этот код работать?
