HPD регион из эмпирического образца? Как отобразить lb и ub?

Question

Ниже приведена функция R, сгенерированная для вычисления интервала HPD 95% из неизвестного унимодального апостериорного распределения с помощью pdf f. (HDP95). Запустив этот код, как мне на самом деле сгенерировать эмпирическую оценку нижней и верхней границ (то есть что-то похожее на proc print в sas ?!)? Спасибо !!

HPD95<-function(LargeSample){  
order<-sort(LargeSample)  
size<-length(LargeSample)  
  isize<-round(0.95*size)  
  lb<-vector("list", (size-isize))  
  ub<-vector("list", (size-isize))  
  inte<-vector("list", (size-isize))  
  for (i in 1:(size-isize)){  
    lb[[i]]<-order[i]  
    ub<-order[i+isize]  
    inte[[i]]<- ub[[i]]-lb[[i]]  
  }  
  minvalue<-min(unlist(inte))  
  position<-which(inte==minvalue)  
  interval<-c(lb[[position]],ub[[position]])  
  return(interval)  
}

Answer 1

Вы можете быть ленивым:

library("coda")
hfun <- function(x,...) coda::HPDinterval(as.mcmc(matrix(x)),...)

(необходим дополнительный механизм, поскольку HPDinterval ожидает объект mcmc; эта функция применяет его к числовому вектору).

Пример:

set.seed(101); hfun(rgamma(1000,shape=2,scale=1))
##           lower    upper
## var1 0.02420174 4.977631
## attr(,"Probability")
## [1] 0.95