Welcome to ClientVPS Mirrors

example

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

Need a high-speed mirror for your open-source project?
Contact our mirror admin team at info@clientvps.com.

This archive is provided as a free public service to the community.
Proudly supported by infrastructure from VPSPulse , RxServers , BuyNumber , UnitVPS , OffshoreName and secure payment technology by ArionPay.