2017-12-28 4 views
0

부트 스트래핑을 사용하여 95 % 신뢰 구간 경계에 감마 분포의 표본 크기 효과를 표시하려고합니다. 이제 4 개의 서로 다른 표본 크기의 결과를 단일 boxplot에 모아야합니다. 는 R 코드는 다음과 같습니다95 % 신뢰도 결과를 비교하기위한 R 박스 플롯 다른 샘플 크기의 간격

y <- rgamma(30,1,1) + rnorm(30,0,0.01) 
y60 <- rgamma(60,1,1) + rnorm(60,0,0.01) 
y100 <- rgamma(100,1,1) + rnorm(100,0,0.01) 
y200 <- rgamma(200,1,1) + rnorm(200,0,0.01) 
minusL <- function(params, data) { 
-sum(log(dgamma(data, params[1], params[2]))) 
} 
fit <- nlm(minusL, c(1,1), data=y) 
fit 
gammamedian<-function(data) { 
fit <- nlm(minusL, c(1,1), data=data) 
qgamma(.5, fit$estimate[1], fit$estimate[2]) 
} 
gammamedian(y) 
gammamedian(y60) 
gammamedian(y100) 
gammamedian(y200) 
gengamma<- function(data, params){ 
rgamma(length(data), params[1], params[2])} 
library(boot) 
pbootresults<-boot(y, gammamedian, R=1000, sim="parametric",    
ran.gen=gengamma, mle=fit$estimate) 
pbootresults 
boot.ci(pbootresults, type=c("basic", "perc", "norm")) 
pbootresults<-boot(y60, gammamedian, R=1000, sim="parametric", 
ran.gen=gengamma, mle=fit$estimate) 
pbootresults 
boot.ci(pbootresults, type=c("basic", "perc", "norm")) 
pbootresults<-boot(y100, gammamedian, R=1000, sim="parametric",  
ran.gen=gengamma, mle=fit$estimate) 
pbootresults 
boot.ci(pbootresults, type=c("basic", "perc", "norm")) 
pbootresults<-boot(y200, gammamedian, R=1000, sim="parametric", 
ran.gen=gengamma, mle=fit$estimate) 
pbootresults 
boot.ci(pbootresults, type=c("basic", "perc", "norm")) 

[An Excel image example ][1] 


[1]: https://i.stack.imgur.com/JXR6P.jpg 
+0

첫 번째 함수'gammamedian'에'minusL'라는 객체가 무엇입니까? –

+0

상기시켜 줘서 고마워. 그냥 빠져 나갔습니다. 함수 minusL은 주어진 파라미터 값에 대해 에 대한 대수 우도 함수를 (빼기로) 계산하고,이를 최소화하기 위해 nlm을 사용하여 최우 추정을 찾습니다. – Reza

+0

귀하의 질문에 대한 답변입니까? 제발 해결해주세요. –

답변

0

당신은 geom_rect()

가 나는 경우이하고 위해 할 수있는 유사 대한 간단한 예를주고 사용 후 데이터 테이블에 minmax 값을 저장해야 휴식을 취하고 아테네 틱스와 함께 더 나은 경치를 즐기십시오.

library(boot) 
pbootresults<-boot(y, gammamedian, R=1000, sim="parametric",    
        ran.gen=gengamma, mle=fit$estimate) 
pbootresults 
temp2 <-boot.ci(pbootresults, type=c("basic", "perc", "norm")) # store teh results 
pbootresults<-boot(y60, gammamedian, R=1000, sim="parametric", 
        ran.gen=gengamma, mle=fit$estimate) 
pbootresults 
boot.ci(pbootresults, type=c("basic", "perc", "norm")) 
pbootresults<-boot(y100, gammamedian, R=1000, sim="parametric",  
        ran.gen=gengamma, mle=fit$estimate) 
pbootresults 
boot.ci(pbootresults, type=c("basic", "perc", "norm")) 
pbootresults<-boot(y200, gammamedian, R=1000, sim="parametric", 
        ran.gen=gengamma, mle=fit$estimate) 
pbootresults 
temp <- boot.ci(pbootresults, type=c("basic", "perc", "norm")) # store the results 

datatable <- temp$percent %>% as.data.table() 
datatable <- rbind(datatable, temp2$percent %>% as.data.table()) 
datatable[, group :=.I] 

ggplot(data = datatable) + geom_rect(mapping = aes(xmin = group - 0.25, xmax = group+ 0.25, ymin = V4, ymax = V5)) 

enter image description here