2010-04-21 4 views
0

나는 최대 우도 방법을 사용하여 R의 역학 모델을 연구하는 학생입니다. 내 부정적인 로그 우도 함수를 만들었습니다. 그것은 찾고 총의 종류,하지만 여기있다 : 걸리는 대부분의 데이터가 입력이 있기 때문에 첫 번째 줄이 너무 끔찍한 보이는R에서 단 하나 헤세 인 행렬을 반전시키지 않고 신뢰 구간을 얻으려면 어떻게해야합니까?

NLLdiff = function(v1, CV1, v2, CV2, st1 = (czI01 - czV01), st2 = (czI02 - czV02), st01 = czI01, st02 = czI02, tt1 = czT01, tt2 = czT02) { 
    prob1 = (1 + v1 * CV1 * tt1)^(-1/CV1) 
    prob2 = (1 + v2 * CV2 * tt2)^(-1/CV2) 
    -(sum(dbinom(st1, st01, prob1, log = T)) + sum(dbinom(st2, st02, prob2, log = T))) 
} 

이유입니다. 예를 들어 czI01이 이미 선언되었습니다. 나는이 함수를 호출하여 함수에 끔찍한 벡터가 생기지 않도록하기 위해 간단하게했다.

그런 다음 mle2 (라이브러리 bbmle)를 사용하여 CV1, CV2, v1 및 v2에 최적화되었습니다. 이는 조금 더보기 좋고 모양은 다음과 같습니다.

ml.cz.diff = mle2 (NLLdiff, start=list(v1 = vguess, CV1 = cguess, v2 = vguess, CV2 = cguess), method="L-BFGS-B", lower = 0.0001) 

이제 여기까지 모든 것이 잘 작동합니다. ml.cz.diff는 내 데이터에 합리적으로 맞는 플롯으로 바뀔 수있는 가치를 제공합니다. 나는 또한 여러 모델을 가지고 있으며 그것들을 비교하기 위해 AICc 값을 얻을 수있다. 그러나 v1, CV1, v2 및 CV2에서 신뢰 구간을 얻으려고하면 문제가 발생합니다. 기본적으로, 나는 CV1에 대해 부정적인 결과를 얻었는데, 이것은 생물학적 모델에서의 제곱 수와 몇 가지 경고를 실제로 나타 내기 때문에 불가능합니다.

신뢰 구간을 얻는 더 좋은 방법이 있습니까? 또는, 정말로, 여기서 신뢰 구간을 얻는 방법?

내가 우연히 볼 때, 헤센 행렬은 최적화 공간의 일부 값에 대해 단수입니다. 그러나 4 개 이상의 변수를 최적화하고 지나치게 광범위한 프로그래밍 지식이 없기 때문에 헤센에 의존하지 않는 최적의 최적화 방법을 생각해 낼 수 없습니다. 나는 그 문제를 봤어. 내 모델은 나쁘다고 제안 했어.하지만 내 모델이 정말로 끔찍하지 않다는 것을 보여주는 몇 가지 작업을 재구성하고있다. (나는 원래 작품의 플롯처럼 ml.cz.diff를 사용하여 만든 플롯이다.). 나는 또한 설명서의 관련 부분과 볼커의 책 생태 모델 R을 읽었습니다. 나는 또한 다른 최적화 방법을 시도했다. 결과적으로 실행 시간은 길지 만 동일한 오류가 발생했다. "SANN"방법은 한 시간 이내에 실행을 마쳤지 않으므로 결과를보기 위해 기다리지 않았습니다.

간단히 말해서 : 내 신뢰 구간이 나쁘다. R에서 그 (것)들을 고치는 상대적으로 간단한 방법 있는가?

v = -log((c(czI01, czI02) - c(czV01, czV02))/c(czI01, czI02))/c(czT01, czT02) 
vguess = mean(v) 
cguess = var(v)/vguess^2 

그것은 내가 다른 완전히 뭔가를 잘못하고 있어요 가능성도 있습니다,하지만 내 결과는 I 피난처 그래서 합리적인 것 :

czT01 = c(5, 5, 5, 5, 5, 5, 5, 25, 25, 25, 25, 25, 25, 25, 50, 50, 50, 50, 50, 50, 50) 
czT02 = c(5, 5, 5, 5, 5, 10, 10, 10, 10, 10, 25, 25, 25, 25, 25, 50, 50, 50, 50, 50, 75, 75, 75, 75, 75) 
czI01 = c(25, 24, 22, 22, 26, 23, 25, 25, 25, 23, 25, 18, 21, 24, 22, 23, 25, 23, 25, 25, 25) 
czI02 = c(13, 16, 5, 18, 16, 13, 17, 22, 13, 15, 15, 22, 12, 12, 13, 13, 11, 19, 21, 13, 21, 18, 16, 15, 11) 
czV01 = c(1, 4, 5, 5, 2, 3, 4, 11, 8, 1, 11, 12, 10, 16, 5, 15, 18, 12, 23, 13, 22) 
czV02 = c(0, 3, 1, 5, 1, 6, 3, 4, 7, 12, 2, 8, 8, 5, 3, 6, 4, 6, 11, 5, 11, 1, 13, 9, 7) 

나는에 의해 나의 추측을 얻을 :

내 벡터는 그것을 잡지 못했다.

답변

6

제약 조건을 항상 만족하도록 매개 변수화를 변경할 수 있습니다. 가능성을 ln (CV1)과 ln (CV2)의 함수로 다시 작성하면 CV1과 CV2가 완전히 양의 값을 유지할 수 있습니다.

NLLdiff_2 = function(v1, lnCV1, v2, lnCV2, st1 = (czI01 - czV01), st2 = (czI02 - czV02), st01 = czI01, st02 = czI02, tt1 = czT01, tt2 = czT02) { 
prob1 = (1 + v1 * exp(lnCV1) * tt1)^(-1/exp(lnCV1)) 
prob2 = (1 + v2 * exp(lnCV2) * tt2)^(-1/exp(lnCV2)) 
-(sum(dbinom(st1, st01, prob1, log = T)) + sum(dbinom(st2, st02, prob2, log = T))) 
} 
관련 문제