[S] SUMMARY: 95 % Hightest posterior density region

Shyh-Forn Guo (guo@calshp.cals.wisc.edu)
Mon, 24 Aug 1998 09:38:01 -0500

Thanks to Kevin Brand,Dave Krantz and Jim Robison-Cox.
I have forwarded all your replies to Ruby Chang.

For those who are interested with the same question, I make
a brief summary and attach their original mails in below.

95% HPD region:
1, If the density is symmetric and has only one peak (like normal):
quantile(posterior,c(0.025,0.975)) can do the job.
2. If the density has a known form like Beta:
Dave provides us a function to do it.
3. For other form unknown distributions ( not symmetric and has
more than one peaks), build a function is necessary. A function
for general case is not available yet and this is what Ruby interested.
There is no option in CODA to do this either.

At 05:21 PM 8/19/98 -0400, you wrote:
>This maybe a bit naive, but a first cut approach for the univariate case
>might be to use the quantile command,
>> quantile(posterior,c(0.025,0.975))

At 07:04 PM 8/19/98 -0400, you wrote:
>I have a fairly crude function (but it works!) that finds highest
>density intervals for beta distributions. It is based on the
>principle that the density must be the same at the upper and lower
>endpoints of a highest density interval. My procedure uses
>the Newton-Raphson method to approximate the zero for the absolute
>value of the logarithm of the ratio of the densities at the two
>endpoints of the interval. I used the lower-tail probability
>of the beta distribution as the parameter, and used a secant slope
>to estimate the derivative of the absolute log density ratio.
>The function was written to take a 2xN matrix of beta parameters,
>representing N betas, and it returns a 2xN matrix of interval
>endpoints. I don't think it would be hard to adapt this to other
>parametric families of densities for which Splus has d- and q-
>functions. [I use dbeta() to get densities from quantiles, and
>I use qbeta() to get quantiles from the current estimate of lower
>tail probability.]
>Dave Krantz

> Ruby,
> I hope you are not pulling a trick on someone. It would be better to
>write to S-news from your own account instead of someone elses.
> What kind of posterior distributions do you wish to use? If you
>conjugate priors, then the builtin cdf functions will help you find the
>HPD regions, (differences of cdf's give probability of intervals) but if
>you use other priors, your posteriors will not be nice
>and you will have to use some sampling from the distributions of interest.
>If my posterior was normal N(m,s^2) I could approximate 95% HPDR this way:
>> m <- 3;s <- 5
>> x <- seq(m - 3*s,m+3*s,length=150)
>> density.x <- dnorm(x,m,s)
>> den.order <- order(density.x)
>> target <- .95
>> cdf.values <- pnorm(x[den.order],m,s)
>> x[den.order][abs(diff(cdf.values)) <= target][1:2]
>[1] -6.7651 12.7651
>##Check to see if difference in cdf (pnorm) is about .95:
>> diff(pnorm(c(-6.7651, 12.7651),m,s))
>[1] 0.9491828
>So the interval (-6.77, 12.77) is about right.
>But this assumes than the region is connected (ie. the posterior is
>unimodal. In general I doubt there is any nice solution, it will have to
>invlove looping.
>Jim Robison-Cox ____________
>Department of Math Sciences | | phone: (406)994-5340
>2-214 Wilson Hall \ BZN, MT | FAX: (406)994-1789
>Montana State University | *_______|
>Bozeman, MT 59717 \_| e-mail: jimrc@math.montana.edu
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news