Skip to contents

Bayesian HPD Regions

Usage

hpdi(object, ...)

# S4 method for numeric
hpdi(object, level = 0.95, CE = TRUE, ...)

# S4 method for MCMC
hpdi(object, level = 0.95, ...)

Arguments

object

A numeric vector or an MCMC object containing the output of the MCMC algorithm for the parameter.

...

Extra arguments to be passed to stats::density().

level

A length-one numeric vector giving the confidence level.

CE

A logical scalar: are the data expressed in BC/AD years?

Value

A list of numericmatrix giving the lower and upper boundaries of the HPD interval and associated probabilities.

References

Hyndman, R. J. (1996). Computing and graphing highest density regions. American Statistician, 50: 120-126. doi:10.2307/2684423 .

See also

stats::density()

Other statistics: credible(), sensitivity(), summary()

Author

A. Philippe, M.-A. Vibet, T. S. Dye, N. Frerebeau

Examples

## Coerce to MCMC
eve <- as_events(events, calendar = "CE", iteration = 1)
eve <- eve[1:10000, ]

## CE
credible(eve, level = 0.95) # Credible interval
#> $E1
#>         lower    upper    p
#> [1,] -1045.21 -202.337 0.95
#> 
#> $E2
#>         lower    upper    p
#> [1,] -1979.81 -1613.59 0.95
#> 
#> $E3
#>         lower   upper    p
#> [1,] -809.057 -459.09 0.95
#> 
#> $E4
#>         lower    upper    p
#> [1,] -1400.04 -1064.29 0.95
#> 
#> attr(,"calendar")
#> [1] "CE"
hpdi(eve, level = 0.95) # HPD interval
#> $E1
#>          lower     upper    p
#> [1,] -1061.227 -193.7855 0.95
#> 
#> $E2
#>          lower     upper    p
#> [1,] -1980.877 -1609.517 0.95
#> 
#> $E3
#>         lower     upper    p
#> [1,] -805.053 -451.8577 0.95
#> 
#> $E4
#>          lower     upper    p
#> [1,] -1406.409 -1066.006 0.95
#> 
#> attr(,"calendar")
#> [1] "CE"

## BP
eve_BP <- CE_to_BP(eve)
credible(eve_BP, level = 0.95) # Credible interval
#> $E1
#>        lower    upper    p
#> [1,] 2994.21 2151.337 0.95
#> 
#> $E2
#>        lower   upper    p
#> [1,] 3928.81 3562.59 0.95
#> 
#> $E3
#>         lower   upper    p
#> [1,] 2758.057 2408.09 0.95
#> 
#> $E4
#>        lower   upper    p
#> [1,] 3349.04 3013.29 0.95
#> 
#> attr(,"calendar")
#> [1] "BP"
hpdi(eve_BP, level = 0.95) # HPD interval
#> $E1
#>         lower    upper    p
#> [1,] 3010.227 2142.786 0.95
#> 
#> $E2
#>         lower    upper    p
#> [1,] 3929.877 3558.517 0.95
#> 
#> $E3
#>         lower    upper    p
#> [1,] 2754.053 2400.858 0.95
#> 
#> $E4
#>         lower    upper    p
#> [1,] 3355.409 3015.006 0.95
#> 
#> attr(,"calendar")
#> [1] "BP"