2016-09-29 7 views
2

나는 (확률 적) 순수 출생 과정에서 출생률을 추론하기 위해 JAGS를 사용하려고합니다. 화학의 언어 JAGS를 이용한 순수 출생 과정 추론

는이 모델의 반응에 상당 : X-> 2X 레이트 알파 *의 X와

이것은 R 부호 I는 (또한 연쇄 반응의 모델로 간주 될 수 있음) '고정 된 시간에 프로세스를 생성하는 데 사용하고 매개 변수 알파의 유추를 만들기 위해 jags 코드를 사용합니다. [말 (람다 나는 새로운 변수에 알파 * y를 [해제 i-1] 퍼팅처럼 다른 일을 시도

Error in jags.model(textConnection(model_string), data = list(y = y, N = N), : 
    RUNTIME ERROR: 
Compilation error on line 4. 
y[2] is a logical node and cannot be observed 

: 나는이 코드를 실행하면

library(rjags) 

y <- 1; # Starting number of "individuals" 
N <- 25 # number of time samplings 
alpha <- 0.2 # per-capita birth rate 
# Generate the time series 
for(i in 2:N) { 
    y<- c(y,y[i-1]+rpois(1,alpha*y[i-1])) 
}; 

# The jags code 
model_string <- "model{ 
    for(i in 2:N) { 
    New[i] ~ dpois(alpha*y[i-1]) 
    y[i] <- y[i-1] + New[i] 
    } 
    alpha ~ dunif(0, 2) 
}" 

# Create and run the jags model 
model <- jags.model(textConnection(model_string), data = list(y = y,N = N), n.chains = 3, n.adapt= 10000) 
update(model, 5000); # Burnin for 10000 samples 
mcmc_samples <- coda.samples(model, variable.names=c("alpha"), n.iter=5000) 

, 나는 다음과 같은 오류가 i]) New [i-1]에 의해 New [i]에 대한 호출을 변경했지만 아무것도 작동하지 않았습니다. 어떤 생각이 왜 실패할까요? 이것을하는 또 다른 똑똑한 방법?

미리 감사드립니다.

+0

나는 그 답을 다른 곳에서 발견했다. https://sourceforge.net/p/mcmc-jags/discussion/610036/thread/6b159634/ –

답변

1

또 다른 솔루션은 데이터를 시뮬레이트하는 방법을 변경하고 모델과 링크 기능을 사용하는 것입니다.

N <- 25 # number of time samplings 
alpha <- 0.2 # log per-capita birth rate 

# Generate the time series 
steps <- 1:N # simulating 25 steps 
log.y<- alpha*steps # the log-scale expected count 
expected.y <- exp(log.y) # back to the real scale 
y <- rpois(N, expected.y) # add Poisson noise to your expected.y 

# The jags code 
model_string <- "model{ 
    for(i in 1:N) { 
    y[i] ~ dpois(lambda[i]) 
    log(lambda[i]) <- log.lambda[i] 
    log.lambda[i] <- alpha * i 
    } 
    alpha ~ dunif(-10, 10) 
}" 

# Create and run the jags model 
model <- jags.model(textConnection(model_string),inits = list(alpha = 1), data = list(y = y,N = N), n.chains = 1, n.adapt= 10000) 
update(model, 5000); # Burnin for 10000 samples 
mcmc_samples <- coda.samples(model, variable.names=c("alpha"), n.iter=5000) 

아래에서 알 수 있듯이이 모델은 매개 변수 (alpha = 0.2)도 올바르게 검색합니다.

enter image description here

그 지수를 촬영하는 것은 당신에게 출산율 (즉 exp(0.2) = 1.22)를 줄 것이다, 또는 당신이 모델 내에서 그것을하고 alpha 단지 기하 급수적 인 파생 매개 변수를 추적 할 수 있습니다. 이 모델은 다음과 같다 :

model_string <- "model{ 
    for(i in 1:N) { 
    y[i] ~ dpois(lambda[i]) 
    log(lambda[i]) <- log.lambda[i] 
    log.lambda[i] <- alpha * i 
    } 
    alpha ~ dunif(-10, 10) 
    birth.rate <- exp(alpha) 
}" 

그리고 당신은 단지 variable.names 인수에 birth.rate을 추적 할 것입니다.

+0

이것은 분명히 내가 예상 한 것보다 나은 해결책입니다. 그러나 나는 왜 비율 람다가 단순히 알파 *의 지수가 아닌지를 완전히 이해하지 못합니다. –

+1

이 모델은 로그 링크를 사용하기 때문에. 로그 링크의 역수는 지수 함수이며 해석하기가 쉬운 비율로 다시 이동합니다 (비율). [포아송 회귀] (http://en.wikipedia.org/wiki/Poisson_regression)의 Wikipedia 페이지에는 이에 대한 충분한 정보가 있습니다. 예를 들어, 1.2의 증가율 (한 ​​단계에서 다음 단계로 20 % 증가)을 시뮬레이트하고 싶다면 로그의 알파 - 로그 (1.2)를 취하여 시뮬레이션에서 사용하십시오. 이 모델은 또한 공변량을 포함 할 수 있지만 다른 솔루션은 그렇지 않습니다. –

+0

지금 삭제하십시오. 다시 한번 감사드립니다. –