Skip to content

Commit 4646f88

Browse files
nfrerebeaunfrerebeau
committed
Add CredibleInterval() as an internal function
1 parent 68aa1eb commit 4646f88

8 files changed

+68
-36
lines changed

‎DESCRIPTION

+1-2
Original file line numberDiff line numberDiff line change
@@ -27,8 +27,7 @@ Imports:
2727
KernSmooth (>= 2.23),
2828
rjags (>= 4-13),
2929
runjags (>= 2.2.1),
30-
Luminescence (>= 0.9.21),
31-
ArchaeoPhases (>= 1.8)
30+
Luminescence (>= 0.9.21)
3231
Suggests:
3332
testthat (>= 3.1.7),
3433
R.rsp (>= 0.45.0)

‎R/AgeC14_Computation.R

+2-2
Original file line numberDiff line numberDiff line change
@@ -335,8 +335,8 @@ AgeC14_Computation <- function(Data_C14Cal,
335335
cat(paste("Sample name", "\t","Bayes estimate"," Credible interval: \n"))
336336
cat(paste(paste("A_",SampleNames[i],sep=""),"\t",(mean(Sample[,i])),'\n'))
337337
cat("\t\t\t\t\t\t lower bound \t upper bound\n")
338-
HPD_95=ArchaeoPhases::CredibleInterval(Sample[,i],0.95,roundingOfValue=roundingOfValue)
339-
HPD_68=ArchaeoPhases::CredibleInterval(Sample[,i],0.68,roundingOfValue=roundingOfValue)
338+
HPD_95=CredibleInterval(Sample[,i],0.95,roundingOfValue=roundingOfValue)
339+
HPD_68=CredibleInterval(Sample[,i],0.68,roundingOfValue=roundingOfValue)
340340
cat("\t\t\t\t at level 95% \t",(c(HPD_95[2])),"\t\t",(c(HPD_95[3])),"\n")
341341
cat("\t\t\t\t at level 68% \t",(c(HPD_68[2])),"\t\t",(c(HPD_68[3])),"\n")
342342
AgePlot95[i,] <- HPD_95

‎R/AgeS_Computation.R

+6-6
Original file line numberDiff line numberDiff line change
@@ -691,9 +691,9 @@ AgeS_Computation <- function(
691691
'\n'
692692
))
693693
cat("\t\t\t\t\t\t lower bound \t upper bound\n")
694-
HPD_95 = ArchaeoPhases::CredibleInterval(sample[, i], 0.95, roundingOfValue =
694+
HPD_95 = CredibleInterval(sample[, i], 0.95, roundingOfValue =
695695
roundingOfValue)
696-
HPD_68 = ArchaeoPhases::CredibleInterval(sample[, i], 0.68, roundingOfValue =
696+
HPD_68 = CredibleInterval(sample[, i], 0.68, roundingOfValue =
697697
roundingOfValue)
698698
cat(
699699
"\t\t\t\t at level 95% \t",
@@ -731,9 +731,9 @@ AgeS_Computation <- function(
731731
'\n'
732732
))
733733
cat("\t\t\t\t\t\t lower bound \t upper bound\n")
734-
HPD_95 = ArchaeoPhases::CredibleInterval(sample[, (Nb_sample + i)], 0.95, roundingOfValue =
734+
HPD_95 = CredibleInterval(sample[, (Nb_sample + i)], 0.95, roundingOfValue =
735735
roundingOfValue)
736-
HPD_68 = ArchaeoPhases::CredibleInterval(sample[, (Nb_sample + i)], 0.68, roundingOfValue =
736+
HPD_68 = CredibleInterval(sample[, (Nb_sample + i)], 0.68, roundingOfValue =
737737
roundingOfValue)
738738
cat(
739739
"\t\t\t\t at level 95% \t",
@@ -768,9 +768,9 @@ AgeS_Computation <- function(
768768
'\n'
769769
))
770770
cat("\t\t\t\t\t\t lower bound \t upper bound\n")
771-
HPD_95 = ArchaeoPhases::CredibleInterval(echantillon[[1]][, (2 * Nb_sample +
771+
HPD_95 = CredibleInterval(echantillon[[1]][, (2 * Nb_sample +
772772
i)], 0.95, roundingOfValue = roundingOfValue)
773-
HPD_68 = ArchaeoPhases::CredibleInterval(echantillon[[1]][, (2 * Nb_sample +
773+
HPD_68 = CredibleInterval(echantillon[[1]][, (2 * Nb_sample +
774774
i)], 0.68, roundingOfValue = roundingOfValue)
775775
cat(
776776
"\t\t\t\t at level 95% \t",

‎R/Age_Computation.R

+6-6
Original file line numberDiff line numberDiff line change
@@ -282,8 +282,8 @@ Age_Computation <- function(
282282
cat(paste("parameter", "\t","Bayes estimate","\t"," Credible interval \n"))
283283
cat("----------------------------------------------\n")
284284
cat(paste("A","\t\t", round(mean(sample[,1]),3),'\n'))
285-
HPD_95=ArchaeoPhases::CredibleInterval(echantillon[[1]][,1],0.95,roundingOfValue = roundingOfValue)
286-
HPD_68=ArchaeoPhases::CredibleInterval(echantillon[[1]][,1],0.68,roundingOfValue = roundingOfValue)
285+
HPD_95=CredibleInterval(echantillon[[1]][,1],0.95,roundingOfValue = roundingOfValue)
286+
HPD_68=CredibleInterval(echantillon[[1]][,1],0.68,roundingOfValue = roundingOfValue)
287287
cat("\t\t\t\t\t\t lower bound \t upper bound\n")
288288
cat("\t\t\t\t at level 95%\t",round(c(HPD_95[2]),2),"\t\t",round(c(HPD_95[3]),roundingOfValue),"\n")
289289
cat("\t\t\t\t at level 68%\t",round(c(HPD_68[2]),2),"\t\t",round(c(HPD_68[3]),roundingOfValue),"\n")
@@ -294,8 +294,8 @@ Age_Computation <- function(
294294
R[1,c(2,4)]=round(HPD_68[2:3],roundingOfValue)
295295

296296
cat(paste("D","\t\t", round(mean(sample[,2]),roundingOfValue),'\n'))
297-
HPD_95=ArchaeoPhases::CredibleInterval(echantillon[[1]][,2],0.95,roundingOfValue = roundingOfValue)
298-
HPD_68=ArchaeoPhases::CredibleInterval(echantillon[[1]][,2],0.68,roundingOfValue = roundingOfValue)
297+
HPD_95=CredibleInterval(echantillon[[1]][,2],0.95,roundingOfValue = roundingOfValue)
298+
HPD_68=CredibleInterval(echantillon[[1]][,2],0.68,roundingOfValue = roundingOfValue)
299299
cat("\t\t\t\t\t\t lower bound \t upper bound\n")
300300
cat("\t\t\t\t at level 95%\t",round(c(HPD_95[2]),roundingOfValue),"\t\t",round(c(HPD_95[3]),roundingOfValue),"\n")
301301
cat("\t\t\t\t at level 68%\t",round(c(HPD_68[2]),roundingOfValue),"\t\t",round(c(HPD_68[3]),roundingOfValue),"\n")
@@ -306,8 +306,8 @@ Age_Computation <- function(
306306

307307
cat("----------------------------------------------\n")
308308
cat(paste("sD","\t\t", round(mean(sample[,3]),3),'\n'))
309-
HPD_95=ArchaeoPhases::CredibleInterval(echantillon[[1]][,3],0.95,roundingOfValue = roundingOfValue)
310-
HPD_68=ArchaeoPhases::CredibleInterval(echantillon[[1]][,3],0.68,roundingOfValue = roundingOfValue)
309+
HPD_95=CredibleInterval(echantillon[[1]][,3],0.95,roundingOfValue = roundingOfValue)
310+
HPD_68=CredibleInterval(echantillon[[1]][,3],0.68,roundingOfValue = roundingOfValue)
311311
cat("\t\t\t\t\t\t lower bound \t upper bound\n")
312312
cat("\t\t\t\t at level 95%\t",round(c(HPD_95[2]),roundingOfValue),"\t\t",round(c(HPD_95[3]),roundingOfValue),"\n")
313313
cat("\t\t\t\t at level 68%\t",round(c(HPD_68[2]),roundingOfValue),"\t\t",round(c(HPD_68[3]),roundingOfValue),"\n")

‎R/Age_OSLC14.R

+2-2
Original file line numberDiff line numberDiff line change
@@ -829,9 +829,9 @@ for (i in 1:results_runjags$args$Nb_sample) {
829829
'\n'
830830
))
831831
cat("\t\t\t\t\t\t lower bound \t upper bound\n")
832-
HPD_95 = ArchaeoPhases::CredibleInterval(Sample[, i], 0.95, roundingOfValue =
832+
HPD_95 = CredibleInterval(Sample[, i], 0.95, roundingOfValue =
833833
roundingOfValue)
834-
HPD_68 = ArchaeoPhases::CredibleInterval(Sample[, i], 0.68, roundingOfValue =
834+
HPD_68 = CredibleInterval(Sample[, i], 0.68, roundingOfValue =
835835
roundingOfValue)
836836
cat(
837837
"\t\t\t\t at level 95% \t",

‎R/Palaeodose_Computation.R

+4-4
Original file line numberDiff line numberDiff line change
@@ -278,8 +278,8 @@ Palaeodose_Computation<-function(
278278
cat(paste("Parameter", "\t","Bayes estimate","\t"," Credible interval \n"))
279279
cat(paste(paste("D_",SampleNames[i],sep=""),"\t",round(mean(sample[,i]),3),'\n'))
280280
cat("\t\t\t\t\t\t lower bound \t upper bound\n")
281-
HPD_95=ArchaeoPhases::CredibleInterval(sample[,i],0.95)
282-
HPD_68=ArchaeoPhases::CredibleInterval(sample[,i],0.68)
281+
HPD_95=CredibleInterval(sample[,i],0.95)
282+
HPD_68=CredibleInterval(sample[,i],0.68)
283283
cat("\t\t\t\t at level 95% \t",round(c(HPD_95[2]),2),"\t\t",round(c(HPD_95[3]),2),"\n")
284284
cat("\t\t\t\t at level 68% \t",round(c(HPD_68[2]),2),"\t\t",round(c(HPD_68[3]),2),"\n")
285285
PalaeodosePlot95[i,]=HPD_95
@@ -293,8 +293,8 @@ Palaeodose_Computation<-function(
293293
cat(paste("\nParameter", "\t","Bayes estimate","\t"," Credible interval \n"))
294294
cat(paste(paste("sD_",SampleNames[i],sep=""),"\t",round(mean(sample[,(Nb_sample+i)]),3),'\n'))
295295
cat("\t\t\t\t\t\t lower bound \t upper bound\n")
296-
HPD_95=ArchaeoPhases::CredibleInterval(sample[,(Nb_sample+i)],0.95)
297-
HPD_68=ArchaeoPhases::CredibleInterval(sample[,(Nb_sample+i)],0.68)
296+
HPD_95=CredibleInterval(sample[,(Nb_sample+i)],0.95)
297+
HPD_68=CredibleInterval(sample[,(Nb_sample+i)],0.68)
298298
cat("\t\t\t\t at level 95% \t",round(c(HPD_95[2]),2),"\t\t",round(c(HPD_95[3]),2),"\n")
299299
cat("\t\t\t\t at level 68% \t",round(c(HPD_68[2]),2),"\t\t",round(c(HPD_68[3]),2),"\n")
300300

‎R/internals.R

+46
Original file line numberDiff line numberDiff line change
@@ -26,3 +26,49 @@
2626
return(l)
2727

2828
}
29+
30+
#' Bayesian Credible Interval
31+
#'
32+
#' Computes the shortest credible interval of the output of the MCMC algorithm
33+
#' for a single parameter
34+
#' @param a_chain Numeric vector containing the output of the MCMC algorithm
35+
#' for the parameter.
36+
#' @param level Probability corresponding to the level of confidence used for
37+
#' the credible interval, default = 0.95.
38+
#' @param roundingOfValue Integer indicating the number of decimal places to be
39+
#' used, default = 0.
40+
#' @details
41+
#' A \eqn{(100 * level)}\% credible interval is an interval that keeps \eqn{N * (1 - level)}
42+
#' elements of the sample outside the interval.
43+
#'
44+
#' The \eqn{(100 * level)}\% credible interval is the shortest of all those intervals.
45+
#' @return
46+
#' A named vector of values containing the confidence level and the endpoints
47+
#' of the shortest credible interval in calendar years (BC/AD).
48+
#' @examples
49+
#' data(Events); attach(Events)
50+
#' CredibleInterval(Event.1)
51+
#' CredibleInterval(Event.12, 0.50)
52+
#' @author A. Philippe, M.-A. Vibet
53+
#' @noRd
54+
CredibleInterval <- function(a_chain, level = 0.95, roundingOfValue = 0) {
55+
sorted_sample <- sort(a_chain) # Ordering the sample
56+
N <- length(a_chain) # Calculation of the sample size
57+
OutSample <- N * (1 - level) # Calculation of the number of data to be outside the interval
58+
59+
# Combinasion of all credible intervals
60+
I = cbind(sorted_sample[1:(OutSample + 1)],
61+
sorted_sample[(N - OutSample):N])
62+
63+
l = I[, 2] - I[, 1] # Length of intervals
64+
i <- which.min(l) # Look for the shortest interval
65+
66+
# Returns the level and the endpoints rounded
67+
return(
68+
c(
69+
"level" = level,
70+
"Credible.Interval.Inf" = round(I[i, 1], digits = roundingOfValue),
71+
"Credible.Interval.Sup" = round(I[i, 2], digits = roundingOfValue)
72+
)
73+
)
74+
}

‎codemeta.json

+1-14
Original file line numberDiff line numberDiff line change
@@ -184,22 +184,9 @@
184184
},
185185
"sameAs": "https://CRAN.R-project.org/package=Luminescence"
186186
},
187-
"10": {
188-
"@type": "SoftwareApplication",
189-
"identifier": "ArchaeoPhases",
190-
"name": "ArchaeoPhases",
191-
"version": ">= 1.8",
192-
"provider": {
193-
"@id": "https://cran.r-project.org",
194-
"@type": "Organization",
195-
"name": "Comprehensive R Archive Network (CRAN)",
196-
"url": "https://cran.r-project.org"
197-
},
198-
"sameAs": "https://CRAN.R-project.org/package=ArchaeoPhases"
199-
},
200187
"SystemRequirements": "JAGS >= 4.3.2 (https://mcmc-jags.sourceforge.io/)"
201188
},
202-
"fileSize": "3463.75KB",
189+
"fileSize": "3465.233KB",
203190
"citation": [
204191
{
205192
"@type": "SoftwareSourceCode",

0 commit comments

Comments
 (0)