Bayesian Credible Interval

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, ...)

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?

Value

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

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
#> $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"