---
title: "example"
author: "Exercised by Chel Hee Lee"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{example}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
bibliography: references.bib
nocite: | 
    @Walley1991, @Bernard2005
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
knitr::opts_chunk$set(fig.width = 6, fig.height=6)
# knitr::opts_chunk$set(tidy.opts=list(width.cutoff=80), tidy=TRUE)
```

```{r setup}
library(imprecise101)
```


## Imprecise Dirichlet Model 

### 4. Bag of Marbles Example 

This example is taken from @Walley1996a. 

#### 4.1. Probabilities of Future Outcomes 

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

```{r}
op <- idm(nj=1, s=2, N=6, k=4)
c(op$p.lower, op$p.upper)
```

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

```{r}
op <- idm(nj=1, s=1, N=6, k=4)
round(c(op$p.lower, op$p.upper),3)
```

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

```{r}
op <- idm(nj=1, s=0, N=6, k=4)
round(c(op$p.lower, op$p.upper),3)
```

Table 1. $P(R|n)$ for different choices of $\Omega$ and $s$ (on page 20)
```{r}
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)
```

For $M=6$ and $s=1$, the CDF are: 

```{r}
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)
```


#### 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$. 

```{r}
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)

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)
```

#### 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. 

```{r}
round(hpd(alpha=3, beta=5, p=0.95),4) # s=2
round(hpd(alpha=3, beta=5, p=0.90),4) # s=2
round(hpd(alpha=3, beta=5, p=0.50),4) # s=2 (required for message of failure)
round(hpd(alpha=2, beta=5, p=0.95),4) # s=1
round(hpd(alpha=2, beta=5, p=0.90),4) # s=1
round(hpd(alpha=2, beta=5, p=0.50),4) # s=1
```

HPD interval, uniform prior $(s=2) [0.0133, 0.5273]$. 
```{r}
x <- pscl::betaHPD(alpha=2, beta=6, p=0.95, plot=FALSE)
round(x,4)
```

#### 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$. 

```{r}
fn <- function(x) choose(7,1)*(1-x)^6
integrate(f=fn, lower=1/2, upper=1)$value

fn <- function(x) dbeta(x, 3, 5)
integrate(f=fn, lower=1/2, upper=1)$value 
```


## Imprecise Beta Model 

### 5. CT vs ECMO Example 

This example is taken from @Walley1996a. 

### 5.5 Deciding when to terminate randomized trials 

```{r out.width="45%"}
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 @Walley1996a.

```{r out.width="45%", echo=FALSE}
tc <- seq(0,1,0.1)
s <- 2

op <- ibm(n=10, m=6, showplot=TRUE, xlab1="z", main1=expression(paste("Fig 1. Posterior ", theta[c], " based on ", s==2)))

op <- ibm(n=9, m=9, showplot=TRUE, xlab1="z", main1=expression(paste("Fig 2. Posterior ", theta[e], " based on ", s==2)))
```


```{r out.width="45%", echo=FALSE}
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=1-cumsum(ymin)/sum(ymin), type="l", ylab="G(z)", xlab="z", main=expression(paste("Fig 3. ", psi, "=", theta[e]-theta[c], " Walley 1996 p.496" )))
points(x=x, y=1-cumsum(ymax)/sum(ymax), type="l")

dmin <- 1-cumsum(ymin/sum(ymin))
dmax <- 1-cumsum(ymax/sum(ymax))
```

Since the imprecision is the area between lower and upper probabilities, $A = \overline{E} - \underline{E}$ =`r sum((dmax - dmin)*0.02)`


## TODO

* Add the examples of decision problem from @Walley1996b. 


## References 
<div id="refs"></div>

## Appendix