Skip to content

Commit a0c19aa

Browse files
audreyyeoCHgithub-actions[bot]dependabot-preview[bot]danielinteractive
authored
120 issues (#142)
Co-authored-by: github-actions <41898282+github-actions[bot]@users.noreply.github.com> Co-authored-by: 27856297+dependabot-preview[bot]@users.noreply.github.com <27856297+dependabot-preview[bot]@users.noreply.github.com> Co-authored-by: Daniel Sabanes Bove <[email protected]>
1 parent b40a03b commit a0c19aa

File tree

10 files changed

+278
-155
lines changed

10 files changed

+278
-155
lines changed

R/dbetabinom.R

Lines changed: 14 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -83,9 +83,20 @@ h_getBetamixPost <- function(x, n, par, weights) {
8383
assert_number(x, lower = 0, upper = n, finite = TRUE)
8484
assert_number(n, lower = 0, finite = TRUE)
8585
assert_matrix(par, min.rows = 1, max.cols = 2, mode = "numeric")
86-
assert_numeric(weights, min.len = 0, len = nrow(par), finite = TRUE)
87-
# We renormalize weights.
88-
weights <- weights / sum(weights)
86+
assert_numeric(
87+
weights,
88+
min.len = 0,
89+
len = nrow(par),
90+
finite = TRUE
91+
)
92+
assert_true(all(weights >= 0))
93+
# Correcting weights that do not sum to 1
94+
if (sum(weights) != 1) {
95+
warning("Weights have been corrected. Advise to review allocated weights")
96+
# We renormalize weights.
97+
weights <- weights / sum(weights)
98+
}
99+
89100
# We now compute updated parameters.
90101
postPar <- par
91102
postPar[, 1] <- postPar[, 1] + x

R/postprob.R

Lines changed: 9 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -68,7 +68,14 @@ postprobBeta <- function(x, n, p, a = 1, b = 1) {
6868
#'
6969
#' @example examples/postprob.R
7070
#' @export
71-
postprob <- function(x, n, p, parE = c(1, 1), weights, betamixPost, log.p = FALSE) {
71+
postprob <- function(
72+
x,
73+
n,
74+
p,
75+
parE = c(1, 1),
76+
weights,
77+
betamixPost,
78+
log.p = FALSE) {
7279
if (missing(betamixPost)) {
7380
assert_flag(log.p)
7481
if (is.vector(parE)) {
@@ -79,7 +86,7 @@ postprob <- function(x, n, p, parE = c(1, 1), weights, betamixPost, log.p = FALS
7986
}
8087
assert_matrix(parE)
8188
if (missing(weights)) {
82-
weights <- rep(1, nrow(parE))
89+
weights <- rep(1 / nrow(parE), nrow(parE))
8390
}
8491
betamixPost <- lapply(
8592
x,

R/predprob.R

Lines changed: 16 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -45,8 +45,7 @@
4545
#'
4646
#' @example examples/predprob.R
4747
#' @export
48-
predprob <- function(x, n, Nmax, p, thetaT, parE = c(1, 1),
49-
weights) {
48+
predprob <- function(x, n, Nmax, p, thetaT, parE = c(1, 1), weights) {
5049
# m = Nmax - n future observations
5150
m <- Nmax - n
5251
if (is.vector(parE)) {
@@ -55,18 +54,29 @@ predprob <- function(x, n, Nmax, p, thetaT, parE = c(1, 1),
5554
parE <- t(parE)
5655
}
5756
if (missing(weights)) {
58-
weights <- rep(1, nrow(parE))
57+
weights <- rep(1 / nrow(parE), nrow(parE))
5958
}
6059
betamixPost <- h_getBetamixPost(x = x, n = n, par = parE, weights = weights)
61-
6260
density <- with(
6361
betamixPost,
6462
dbetabinomMix(x = 0:m, m = m, par = par, weights = weights)
6563
)
66-
assert_numeric(density, lower = 0, upper = 1 + .Machine$double.eps, finite = TRUE, any.missing = FALSE)
64+
assert_numeric(
65+
density,
66+
lower = 0,
67+
upper = 1 + .Machine$double.eps,
68+
finite = TRUE,
69+
any.missing = FALSE
70+
)
6771
assert_number(thetaT, lower = 0, upper = 1, finite = TRUE)
6872
# posterior probabilities to be above threshold p
69-
posterior <- postprob(x = x + c(0:m), n = Nmax, p = p, parE = parE, weights = weights)
73+
posterior <- postprob(
74+
x = x + c(0:m),
75+
n = Nmax,
76+
p = p,
77+
parE = parE,
78+
weights = weights
79+
)
7080
list(
7181
result = sum(density * (posterior > thetaT)),
7282
table = data.frame(

examples/postprob.R

Lines changed: 14 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -10,21 +10,23 @@ postprob(x = 16, n = 23, p = 0.60, par = c(0.6, 0.4))
1010
# 2 component beta mixture prior :
1111
# i.e., P_E ~ 0.6*beta(0.6,0.4) + 0.4*beta(1,1) and Pr(P_E > p | data) = 0.823
1212
postprob(
13-
x = 16, n = 23, p = 0.60,
14-
par =
15-
rbind(
16-
c(0.6, 0.4),
17-
c(1, 1)
18-
),
13+
x = 16,
14+
n = 23,
15+
p = 0.60,
16+
par = rbind(
17+
c(0.6, 0.4),
18+
c(1, 1)
19+
),
1920
weights = c(0.6, 0.4)
2021
)
2122

2223
postprob(
23-
x = 0:23, n = 23, p = 0.60,
24-
par =
25-
rbind(
26-
c(0.6, 0.4),
27-
c(1, 1)
28-
),
24+
x = 0:23,
25+
n = 23,
26+
p = 0.60,
27+
par = rbind(
28+
c(0.6, 0.4),
29+
c(1, 1)
30+
),
2931
weights = c(0.6, 0.4)
3032
)

examples/predprob.R

Lines changed: 21 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -2,25 +2,41 @@
22
# Nmax = 40, x = 16, n = 23, beta(0.6,0.4) prior distribution,
33
# thetaT = 0.9. The control response rate is 60%:
44
predprob(
5-
x = 16, n = 23, Nmax = 40, p = 0.6, thetaT = 0.9,
5+
x = 16,
6+
n = 23,
7+
Nmax = 40,
8+
p = 0.6,
9+
thetaT = 0.9,
610
parE = c(0.6, 0.4)
711
)
812

913
# Lowering/Increasing the probability threshold thetaT of course increases
1014
# /decreases the predictive probability of success, respectively:
1115
predprob(
12-
x = 16, n = 23, Nmax = 40, p = 0.6, thetaT = 0.8,
16+
x = 16,
17+
n = 23,
18+
Nmax = 40,
19+
p = 0.6,
20+
thetaT = 0.8,
1321
parE = c(0.6, 0.4)
1422
)
1523

1624
predprob(
17-
x = 16, n = 23, Nmax = 40, p = 0.6, thetaT = 0.95,
25+
x = 16,
26+
n = 23,
27+
Nmax = 40,
28+
p = 0.6,
29+
thetaT = 0.95,
1830
parE = c(0.6, 0.4)
1931
)
2032

2133
# Mixed beta prior
2234
predprob(
23-
x = 20, n = 23, Nmax = 40, p = 0.6, thetaT = 0.9,
35+
x = 20,
36+
n = 23,
37+
Nmax = 40,
38+
p = 0.6,
39+
thetaT = 0.9,
2440
parE = rbind(c(1, 1), c(25, 15)),
25-
weights = c(3, 1)
41+
weights = c(0.5, 0.5)
2642
)

man/postprob.Rd

Lines changed: 14 additions & 12 deletions
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

man/predprob.Rd

Lines changed: 21 additions & 5 deletions
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

0 commit comments

Comments
 (0)