Skip to content

Commit cf89501

Browse files
FBartosFBartos
authored
Fix release (#218)
* part1 * updates 2 * updates 3 * updates 4 * fixes: jasp-stats/jasp-test-release#2649 * fixes: jasp-stats/jasp-test-release#2651 * fixes: jasp-stats/jasp-test-release#2659 * fixes: jasp-stats/jasp-test-release#2661 * as before for estimation * fixes: jasp-stats/jasp-test-release#2663 * fixes
1 parent f00d000 commit cf89501

11 files changed

+332
-148
lines changed

‎R/LSbinomialcommon.R

Lines changed: 52 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -216,13 +216,15 @@
216216

217217
output <- list(
218218
distribution = distributionText,
219-
mean = .betaMeanLS(alpha, beta, lower, upper),
220-
median = .qbetaLS(.5, alpha, beta, lower, upper),
221-
mode = .betaModeLS(alpha, beta, lower, upper),
222-
lCI = .qbetaLS(.025, alpha, beta, lower, upper),
223-
uCI = .qbetaLS(.975, alpha, beta, lower, upper)
219+
mean = try(.betaMeanLS(alpha, beta, lower, upper)),
220+
median = try(.qbetaLS(.5, alpha, beta, lower, upper)),
221+
mode = try(.betaModeLS(alpha, beta, lower, upper)),
222+
lCI = try(.qbetaLS(.025, alpha, beta, lower, upper)),
223+
uCI = try(.qbetaLS(.975, alpha, beta, lower, upper))
224224
)
225225

226+
# check and recover from errors
227+
output <- .checkAndRecoverErrors(output)
226228

227229
return(output)
228230
}
@@ -262,11 +264,12 @@
262264
tempPostUper <- stats::pbeta(tempPrior[["truncationUpper"]], data$nSuccesses + tempPrior[["betaPriorAlpha"]], data$nFailures + tempPrior[["betaPriorBeta"]])
263265

264266
logLik[i] <- logLik[i] + log(tempPostUper - tempPostLower) - log(tempPriorUpper - tempPriorLower)
265-
266267
}
267268

268269
}
269270

271+
if (is.infinite(logLik[i]))
272+
logLik[i] <- NA
270273

271274
} else {
272275
logLik[i] <- 0
@@ -277,7 +280,7 @@
277280
if (data$nSuccesses + data$nFailures > 0) {
278281

279282
priorWeightLogLik <- log(prior) + logLik
280-
normConst <- log(sum(exp(priorWeightLogLik)))
283+
normConst <- log(sum(exp(priorWeightLogLik), na.rm = TRUE))
281284
posterior <- exp(priorWeightLogLik - normConst)
282285

283286
} else {
@@ -335,14 +338,17 @@
335338

336339
output <- list(
337340
distribution = distributionText,
338-
mean = .bbinomMeanLS(n, alpha, beta, lower, upper) / d,
339-
median = .qbbinomLS(.5, n, alpha, beta, lower, upper) / d,
340-
mode = .bbinomModeLS(n, alpha, beta, lower, upper, prop = prop),
341-
lCI = .qbbinomLS(0.025, n, alpha, beta, lower, upper) / d,
342-
uCI = .qbbinomLS(0.975, n, alpha, beta, lower, upper) / d,
343-
SD = .bbinomSdLS(n, alpha, beta, lower, upper) / d
341+
mean = try(.bbinomMeanLS(n, alpha, beta, lower, upper) / d),
342+
median = try(.qbbinomLS(.5, n, alpha, beta, lower, upper) / d),
343+
mode = try(.bbinomModeLS(n, alpha, beta, lower, upper, prop = prop)),
344+
lCI = try(.qbbinomLS(0.025, n, alpha, beta, lower, upper) / d),
345+
uCI = try(.qbbinomLS(0.975, n, alpha, beta, lower, upper) / d),
346+
SD = try(.bbinomSdLS(n, alpha, beta, lower, upper) / d)
344347
)
345348

349+
# check and recover from errors
350+
output <- .checkAndRecoverErrors(output)
351+
346352
return(output)
347353
}
348354
}
@@ -666,10 +672,11 @@
666672
if (lower == 0 && upper == 1) {
667673
return(stats::dbeta(x, alpha, beta))
668674
} else {
669-
num <- stats::dbeta(x, alpha, beta)
670-
den <- pbeta(upper, alpha, beta) - pbeta(lower, alpha, beta)
671-
lik <- num/den
675+
num <- stats::dbeta(x, alpha, beta, log = TRUE)
676+
den <- log(pbeta(upper, alpha, beta) - pbeta(lower, alpha, beta))
677+
lik <- exp(num - den)
672678
lik[x < lower | x > upper] <- 0
679+
lik[is.nan(lik)] <- NA
673680
return(lik)
674681
}
675682
}
@@ -683,6 +690,7 @@
683690
p <- (p - C1) / (C2 - C1)
684691
p[x < lower] <- 0
685692
p[x > upper] <- 1
693+
p[is.nan(p)] <- NA
686694
return(p)
687695
}
688696
}
@@ -692,26 +700,29 @@
692700
} else {
693701

694702
q <- sapply(x, function(xi) {
695-
stats::optim(
703+
tempOptim <- stats::optim(
696704
par = (lower + upper) / 2 ,
697705
fn = function(p) ( .pbetaLS(p, alpha, beta, lower, upper) - xi)^2,
698706
lower = lower,
699707
upper = upper,
700-
method = "L-BFGS-B",
701-
control = list(
702-
factr = 1e3
703-
)
704-
)$par
708+
method = "Brent"
709+
)
710+
if (is.na(tempOptim$value) || tempOptim$value > 1e-5)
711+
return(NA)
712+
else
713+
return(tempOptim$par)
705714
})
706-
707715
return(q)
708716
}
709717
}
710718
.betaMeanLS <- function(alpha, beta, lower = 0, upper = 1){
711719
if (lower == 0 && upper == 1) {
712720
return(alpha / (alpha + beta))
713721
} else {
714-
return(stats::integrate(function(x) x * .dbetaLS(x, alpha, beta, lower, upper), lower, upper)$value)
722+
tempMean <- stats::integrate(function(x) x * .dbetaLS(x, alpha, beta, lower, upper), lower, upper)$value
723+
if (tempMean < lower || tempMean > upper)
724+
stop("Mean is outside of the truncation bounds -- numerical precision error.")
725+
return(tempMean)
715726
}
716727
}
717728
.betaModeLS <- function(alpha, beta, lower = 0, upper = 1){
@@ -746,12 +757,16 @@
746757
return(extraDistr::dbbinom(x, size, alpha, beta))
747758
} else {
748759
return(sapply(x, function(xi) {
749-
integrate(
760+
tempInt <- try(integrate(
750761
f = function(p) {
751762
stats::dbinom(x = xi, size = size, prob = p) * .dbetaLS(p, alpha, beta, lower, upper)
752763
},
753764
lower = lower,
754-
upper = upper)$value
765+
upper = upper))
766+
if (jaspBase::isTryError(tempInt))
767+
return(NA)
768+
else
769+
return(tempInt$value)
755770
}))
756771
}
757772
}
@@ -892,8 +907,8 @@
892907
x <- qbinom(c((1 - coverage)/2 + 1e-5, 1 - (1 - coverage)/2), n, prior[["spikePoint"]])
893908
else if (prior[["type"]] == "beta")
894909
x <- c(
895-
.qbbinomLS((1 - coverage)/2 + 1e-5, n, prior[["betaPriorAlpha"]] + data$nSuccesses, prior[["betaPriorBeta"]] + data$nFailures, prior[["truncationLower"]], prior[["truncationUpper"]]),
896-
.qbbinomLS(1 - (1 - coverage)/2, n , prior[["betaPriorAlpha"]] + data$nSuccesses, prior[["betaPriorBeta"]] + data$nFailures, prior[["truncationLower"]], prior[["truncationUpper"]])
910+
.qbbinomLS((1 - coverage)/2 + 1e-5, n, prior[["betaPriorAlpha"]] + data$nSuccesses, prior[["betaPriorBeta"]] + data$nFailures, prior[["truncationLower"]], prior[["truncationUpper"]]),
911+
.qbbinomLS(1 - (1 - coverage)/2, n, prior[["betaPriorAlpha"]] + data$nSuccesses, prior[["betaPriorBeta"]] + data$nFailures, prior[["truncationLower"]], prior[["truncationUpper"]])
897912
)
898913

899914

@@ -908,9 +923,14 @@
908923

909924
if (prior[["type"]] == "spike")
910925
coverage <- ifelse (lCI <= prior[["spikePoint"]] & prior[["spikePoint"]] <= uCI, 1, 0)
911-
else if (prior[["type"]] == "beta")
912-
coverage <- .pbetaLS(uCI, prior[["betaPriorAlpha"]] + data$nSuccesses, prior[["betaPriorBeta"]] + data$nFailures, prior[["truncationLower"]], prior[["truncationUpper"]]) -
913-
.pbetaLS(lCI, prior[["betaPriorAlpha"]] + data$nSuccesses, prior[["betaPriorBeta"]] + data$nFailures, prior[["truncationLower"]], prior[["truncationUpper"]])
926+
else if (prior[["type"]] == "beta") {
927+
if (uCI >= prior[["truncationUpper"]] && lCI <= prior[["truncationLower"]])
928+
coverage <- 1
929+
else
930+
coverage <- .pbetaLS(uCI, prior[["betaPriorAlpha"]] + data$nSuccesses, prior[["betaPriorBeta"]] + data$nFailures, prior[["truncationLower"]], prior[["truncationUpper"]]) -
931+
.pbetaLS(lCI, prior[["betaPriorAlpha"]] + data$nSuccesses, prior[["betaPriorBeta"]] + data$nFailures, prior[["truncationLower"]], prior[["truncationUpper"]])
932+
}
933+
914934

915935

916936
} else if (type == "prediction") {
@@ -924,6 +944,7 @@
924944

925945
}
926946

947+
927948
dat <- data.frame(xStart = lCI, xEnd = uCI, g = "custom", coverage = coverage, parameter = "theta")
928949
return(dat)
929950
}

0 commit comments

Comments
 (0)