library(imprecise101)
This example is taken from Peter Walley (1996).
For s=2, ¯P(R|n)=0.375 and P_(R|n)=0.125.
<- idm(nj=1, s=2, N=6, k=4)
op c(op$p.lower, op$p.upper)
#> [1] 0.125 0.375
For s=1, ¯P(R|n)=0.286 and P_(R|n)=0.143.
<- idm(nj=1, s=1, N=6, k=4)
op round(c(op$p.lower, op$p.upper),3)
#> [1] 0.143 0.286
For s=0, P(R|n)=0.167 irrespective of Ω.
<- idm(nj=1, s=0, N=6, k=4)
op round(c(op$p.lower, op$p.upper),3)
#> [1] 0.167 0.167
Table 1. P(R|n) for different choices of Ω and s (on page 20)
<- c(idm(nj=1, s=1, N=6, k=4)$p,
r1 idm(nj=1, s=2, N=6, k=4)$p,
idm(nj=1, s=4/2, N=6, k=4)$p,
idm(nj=1, s=4, N=6, k=4)$p) # Omega 1
<- c(idm(nj=1, s=1, N=6, k=2)$p,
r2 idm(nj=1, s=2, N=6, k=2)$p,
idm(nj=1, s=2/2, N=6, k=2)$p,
idm(nj=1, s=2, N=6, k=2)$p) # Omega 2
<- c(idm(nj=1, s=1, N=6, k=3, cA=2)$p,
r3 idm(nj=1, s=2, N=6, k=3, cA=2)$p,
idm(nj=1, s=3/2, N=6, k=3, cA=2)$p,
idm(nj=1, s=3, N=6, k=3, cA=2)$p) # Omega 3
<- rbind(r1, r2, r3)
tb1 rownames(tb1) <- c("Omega1", "Omega2", "Omega3")
colnames(tb1) <- c("s=1", "s=2", "s=k/2", "s=k")
round(tb1,3)
#> s=1 s=2 s=k/2 s=k
#> Omega1 0.179 0.188 0.188 0.200
#> Omega2 0.214 0.250 0.214 0.250
#> Omega3 0.238 0.292 0.267 0.333
For M=6 and s=1, the CDF are:
<- cbind(
mat unlist(pbetabinom(M=6, x=1, s=1, N=6, y=0)),
unlist(pbetabinom(M=6, x=1, s=1, N=6, y=1)),
unlist(pbetabinom(M=6, x=1, s=1, N=6, y=2)),
unlist(pbetabinom(M=6, x=1, s=1, N=6, y=3)),
unlist(pbetabinom(M=6, x=1, s=1, N=6, y=4)),
unlist(pbetabinom(M=6, x=1, s=1, N=6, y=5)),
unlist(pbetabinom(M=6, x=1, s=1, N=6, y=6))
)colnames(mat) <- c("y=0", "y=1", "y=2", "y=3", "y=4", "y=5", "y=6")
round(mat, 3)
#> y=0 y=1 y=2 y=3 y=4 y=5 y=6
#> p.l 0.227 0.500 0.727 0.879 0.960 0.992 1
#> p.u 0.500 0.773 0.909 0.970 0.992 0.999 1
For s=2, ¯E(θR|n)=0.375, E_(θR|n)=0.125, ¯σ(θR|n)=0.188, and σ_(θR|n)=0.110.
For s=1, ¯E(θR|n)=0.286, E_(θR|n)=0.143, ¯σ(θR|n)=0.164, and σ_(θR|n)=0.124.
<- idm(nj=1, s=2, N=6, k=4)
op round(c(op$p.upper, op$p.lower, op$s.upper, op$s.lower),3)
#> [1] 0.375 0.125 0.188 0.110
<- idm(nj=1, s=1, N=6, k=4)
op round(c(op$p.upper, op$p.lower, op$s.upper, op$s.lower),3)
#> [1] 0.286 0.143 0.164 0.124
For s=2, 95%, 90%, and 50% credible intervals are [0.0031,0.6587], [0.0066,0.5962], and [0.0481,0.3656], respectively.
For s=1, 95%, 90%, and 50% credible intervals are [0.0076,0.5834], [0.0150,0.5141], [0.0761,0.2958], respectively.
round(hpd(alpha=3, beta=5, p=0.95),4) # s=2
#> a b
#> 0.0031 0.6588
round(hpd(alpha=3, beta=5, p=0.90),4) # s=2
#> a b
#> 0.0066 0.5962
round(hpd(alpha=3, beta=5, p=0.50),4) # s=2 (required for message of failure)
#> a b
#> 0.3641 0.9048
round(hpd(alpha=2, beta=5, p=0.95),4) # s=1
#> a b
#> 0.0076 0.5834
round(hpd(alpha=2, beta=5, p=0.90),4) # s=1
#> a b
#> 0.015 0.514
round(hpd(alpha=2, beta=5, p=0.50),4) # s=1
#> a b
#> 0.0760 0.2957
HPD interval, uniform prior (s=2)[0.0133,0.5273].
<- pscl::betaHPD(alpha=2, beta=6, p=0.95, plot=FALSE)
x round(x,4)
#> [1] 0.0133 0.5273
Test H0:θR≥1/2 against Ha:θR<1/2. The posterior lower and upper probabilities of H0 are P_(H0|n)=0.00781 and ¯P(H0|n)=0.227.
<- function(x) choose(7,1)*(1-x)^6
fn integrate(f=fn, lower=1/2, upper=1)$value
#> [1] 0.0078125
<- function(x) dbeta(x, 3, 5)
fn integrate(f=fn, lower=1/2, upper=1)$value
#> [1] 0.2265625
This example is taken from Peter Walley (1996).
<- seq(-0.99, 0.99, 0.02)
x <- ymin <- numeric(length(x))
ymax for(i in 1:length(x)) ymin[i] <- dbetadif(x=x[i], a1=9,b1=2,a2=8,b2=4)
for(i in 1:length(x)) ymax[i] <- dbetadif(x=x[i], a1=11,b1=0.01,a2=6,b2=6)
plot(x=x, y=cumsum(ymin)/sum(ymin), type="l", ylab="F(z)", xlab="z",
main=expression(paste("Fig 1. Posterior upper and lower CDFs for ", psi,
"=", theta[e]-theta[c])))
points(x=x, y=cumsum(ymax)/sum(ymax), type="l")
This example is taken from Peter Walley (1996).
Since the imprecision is the area between lower and upper probabilities, A=¯E−E_ =0.3475766