Skip to contents

Computes the shortest credible interval of the output of the MCMC algorithm for a single parameter.

Usage

credible(object, ...)

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

# S4 method for MCMC
credible(object, level = 0.95, simplify = TRUE)

Arguments

object

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

...

Currently not used.

level

A length-one numeric vector giving the confidence level.

CE

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

simplify

A logical scalar: should the results should be simplified?

Value

If simplify is TRUE (the default), returns a data.frame giving the lower and upper boundaries of the credible interval and associated probabilities. Else, returns a list of numeric

matrix.

Details

A \((100 \times level)\) % credible interval is an interval that keeps \(N \times (1 - level)\) elements of the sample outside the interval.

The \((100 \times level)\) % credible interval is the shortest of all those intervals.

For instance, the 95% credible interval is the central portion of the posterior distribution that contains 95% of the values.

See also

Other statistics: hpdi(), 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
#>   name     lower     upper    p
#> 1   E1 -1045.210  -202.337 0.95
#> 2   E2 -1979.810 -1613.590 0.95
#> 3   E3  -809.057  -459.090 0.95
#> 4   E4 -1400.040 -1064.290 0.95
hpdi(eve, level = 0.68) # HPD interval
#>   name      lower      upper    p
#> 1   E1 -1009.3026  -774.1160 0.35
#> 2   E1  -490.0595  -257.9273 0.33
#> 3   E2 -1891.7506 -1690.3908 0.68
#> 4   E3  -757.2008  -604.5292 0.68
#> 5   E4 -1311.7275 -1149.4158 0.68

## BP
eve_BP <- CE_to_BP(eve)
credible(eve_BP, level = 0.95) # Credible interval
#>   name    lower    upper    p
#> 1   E1 2994.210 2151.337 0.95
#> 2   E2 3928.810 3562.590 0.95
#> 3   E3 2758.057 2408.090 0.95
#> 4   E4 3349.040 3013.290 0.95
hpdi(eve_BP, level = 0.95) # HPD interval
#>   name    lower    upper    p
#> 1   E1 3010.227 2142.786 0.95
#> 2   E2 3929.877 3558.517 0.95
#> 3   E3 2754.053 2400.858 0.95
#> 4   E4 3355.409 3015.006 0.95