The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

example

Exercised by Chel Hee Lee

library(imprecise101)

Imprecise Dirichlet Model

4. Bag of Marbles Example

This example is taken from Peter Walley (1996).

4.1. Probabilities of Future Outcomes

For \(s=2\), \(\overline{P}(R|n)=0.375\) and \(\underline{P}(R|n)=0.125\).

op <- idm(nj=1, s=2, N=6, k=4)
c(op$p.lower, op$p.upper)
#> [1] 0.125 0.375

For \(s=1\), \(\overline{P}(R|n)=0.286\) and \(\underline{P}(R|n)=0.143\).

op <- idm(nj=1, s=1, N=6, k=4)
round(c(op$p.lower, op$p.upper),3)
#> [1] 0.143 0.286

For \(s=0\), \(P(R|n)=0.167\) irrespective of \(\Omega\).

op <- idm(nj=1, s=0, N=6, k=4)
round(c(op$p.lower, op$p.upper),3)
#> [1] 0.167 0.167

Table 1. \(P(R|n)\) for different choices of \(\Omega\) and \(s\) (on page 20)

r1 <- c(idm(nj=1, s=1, N=6, k=4)$p, 
        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

r2 <- c(idm(nj=1, s=1, N=6, k=2)$p, 
        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 

r3 <- c(idm(nj=1, s=1, N=6, k=3, cA=2)$p, 
        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

tb1 <- rbind(r1, r2, r3)
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:

mat <- cbind(
  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

4.2 Means and Standard Deviations for \(\theta_R\)

For \(s=2\), \(\overline{E}(\theta_R|n) = 0.375\), \(\underline{E}(\theta_R|n)=0.125\), \(\overline{\sigma}(\theta_R|n)=0.188\), and \(\underline{\sigma}(\theta_R|n)=0.110\).

For \(s=1\), \(\overline{E}(\theta_R|n) = 0.286\), \(\underline{E}(\theta_R|n)=0.143\), \(\overline{\sigma}(\theta_R|n)=0.164\), and \(\underline{\sigma}(\theta_R|n)=0.124\).

op <- idm(nj=1, s=2, N=6, k=4)
round(c(op$p.upper, op$p.lower, op$s.upper, op$s.lower),3)
#> [1] 0.375 0.125 0.188 0.110

op <- idm(nj=1, s=1, N=6, k=4)
round(c(op$p.upper, op$p.lower, op$s.upper, op$s.lower),3)
#> [1] 0.286 0.143 0.164 0.124

4.3 Credible Intervals for \(\theta_R\)

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]\).

x <- pscl::betaHPD(alpha=2, beta=6, p=0.95, plot=FALSE)
round(x,4)
#> [1] 0.0133 0.5273

4.4 Testing Hypotheses about \(\theta_R\)

Test \(H_0: \theta_R \ge 1/2\) against \(H_a: \theta_R < 1/2\). The posterior lower and upper probabilities of \(H_0\) are \(\underline{P}(H_0|n)=0.00781\) and \(\overline{P}(H_0|n)=0.227\).

fn <- function(x) choose(7,1)*(1-x)^6
integrate(f=fn, lower=1/2, upper=1)$value
#> [1] 0.0078125

fn <- function(x) dbeta(x, 3, 5)
integrate(f=fn, lower=1/2, upper=1)$value 
#> [1] 0.2265625

Imprecise Beta Model

5. CT vs ECMO Example

This example is taken from Peter Walley (1996).

5.5 Deciding when to terminate randomized trials

x <- seq(-0.99, 0.99, 0.02)
ymax <- ymin <- numeric(length(x))
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")

Further inferences concerning \(\theta_c\), \(\theta_e\), and \(\psi\)

This example is taken from Peter Walley (1996).

Since the imprecision is the area between lower and upper probabilities, \(A = \overline{E} - \underline{E}\) =0.3475766

TODO

References

Bernard, Jean-Marc. 2005. “An Introduction to the Imprecise Dirichlet Model for Multinomial Data.” International Journal of Approximate Reasoning 39 (2): 123–50. https://doi.org/https://doi.org/10.1016/j.ijar.2004.10.002.
Walley, P. 1991. Statistical Reasoning with Imprecise Probabilities. Chapman & Hall/CRC Monographs on Statistics & Applied Probability. Taylor & Francis.
Walley, Peter. 1996. “Inferences from Multinomial Data: Learning about a Bag of Marbles.” Journal of the Royal Statistical Society. Series B (Methodological) 58 (1): 3–57. http://www.jstor.org/stable/2346164.
Walley, Peter, Lyle Gurrin, and Paul Burton. 1996. “Analysis of Clinical Data Using Imprecise Prior Probabilities.” Journal of the Royal Statistical Society. Series D (The Statistician) 45 (4): 457–85. http://www.jstor.org/stable/2988546.

Appendix

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.