Chapter 3 - Sampling the Imaginary

This chapter introduced the basic procedures for manipulating posterior distributions. Our fundamental tool is samples of parameter values drawn from the posterior distribution. These samples can be used to produce intervals, point estimates, posterior predictive checks, as well as other kinds of simulations. Posterior predictive checks combine uncertainty about parameters, as described by the posterior distribution, with uncertainty about outcomes, as described by the assumed likelihood function. These checks are useful for verifying that your software worked correctly. They are also useful for prospecting for ways in which your models are inadequate.

Place each answer inside the code chunk (grey box). The code chunks should contain a text response or a code that completes/answers the question or activity requested. Make sure to include plots if the question requests them.

Finally, upon completion, name your final output .html file as: YourName_ANLY505-Year-Semester.html and publish the assignment to your R Pubs account and submit the link to Canvas. Each question is worth 5 points

Questions

3-1. Suppose the globe tossing data had turned out to be 8 water in 15 tosses. Construct the posterior distribution, using grid approximation. Plot the posterior. Use the same flat prior as before.

# Determination of Posterior distribution probability with the utilization of grid approximation approach.

Pstr_grid <- seq(from = 0, to = 1, length.out = 1000 )

Prior <- rep(1, 1000)

# Therefore the likelihood of the distribution under the assumption of 8 water in 15 tosses would be as under:

Likelihood <- dbinom(8, size = 15, prob = Pstr_grid)

Posterior <- Likelihood * Prior

Posterior <- Posterior / sum(Posterior)

# Plot construction for the posterior distribution.

Plot_Pstr_Distrbtn<- plot( x = Pstr_grid, y = Posterior, type = "l" )

Plot_Pstr_Distrbtn
## NULL

3-2. Draw 10,000 samples from the grid approximation from above. Then use the samples to calculate the 90% HPDI for p.

library("rstan")

library("rethinking")


# Determination and drawing 10,000 samples from the earlier calculated grid approximation:

Samples <- sample( Pstr_grid, prob = Posterior, size = 10000, replace = TRUE )


# Utilization of the sample for the calculation of the 90 % HPDI for P.

HPDI_Calcultn_90_Prcnt <- HPDI(Samples, prob = 0.90 )


HPDI_Calcultn_90_Prcnt
##      |0.9      0.9| 
## 0.3393393 0.7217217
## Hence 90 % HPDI for p is within [0.33 , 0.72] 

3-3. Construct a posterior predictive check for this model and data. This means simulate the distribution of samples, averaging over the posterior uncertainty in p. What is the probability of observing 8 water in 15 tosses?

# Determination of Posterior Predictive check for the above model and data.

Data <- rbinom(10000, size = 15, prob = Samples)

# Plotting the data histogram

Histgrm_Plot <- hist(Data)

Histgrm_Plot
## $breaks
##  [1]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
## 
## $counts
##  [1]   39  104  270  551  814 1142 1313 1499 1402 1152  876  484  250   86   18
## 
## $density
##  [1] 0.0039 0.0104 0.0270 0.0551 0.0814 0.1142 0.1313 0.1499 0.1402 0.1152
## [11] 0.0876 0.0484 0.0250 0.0086 0.0018
## 
## $mids
##  [1]  0.5  1.5  2.5  3.5  4.5  5.5  6.5  7.5  8.5  9.5 10.5 11.5 12.5 13.5 14.5
## 
## $xname
## [1] "Data"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
# Determination for the probability of 8 water in 15 tosses.

Problty <- length(Data[Data == 8]) / length(Data)

Problty
## [1] 0.1499
# Hence the he probability of 8 water in 15 tosses would be 14.5 %.

3-4. Using the posterior distribution constructed from the new (8/15) data, now calculate the probability of observing 6 water in 9 tosses.

# Determination of the probability of observing 6 water in 9 tosses.

Data_Rvsd <- rbinom(10000, size = 9, prob = Samples)

# Plotting the new data histogram

Histgrm_Plot_Rvsd <- hist(Data_Rvsd)

Histgrm_Plot_Rvsd
## $breaks
##  [1] 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0
## 
## $counts
##  [1]   43  291    0  811    0 1376    0 1918    0 2031    0 1807    0 1144    0
## [16]  475    0  104
## 
## $density
##  [1] 0.0086 0.0582 0.0000 0.1622 0.0000 0.2752 0.0000 0.3836 0.0000 0.4062
## [11] 0.0000 0.3614 0.0000 0.2288 0.0000 0.0950 0.0000 0.0208
## 
## $mids
##  [1] 0.25 0.75 1.25 1.75 2.25 2.75 3.25 3.75 4.25 4.75 5.25 5.75 6.25 6.75 7.25
## [16] 7.75 8.25 8.75
## 
## $xname
## [1] "Data_Rvsd"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
mean (Data_Rvsd == 6)
## [1] 0.1807
## Probability calculation of observing 6 water in 9 tosses


Problty_Rvsd <- length(Data_Rvsd[Data_Rvsd == 6]) / length(Data_Rvsd)

Problty_Rvsd
## [1] 0.1807
# Hence the Probability of observing 6 water in 9 tosses would be 17.2 percent.

3-5. Start over at 3-1, but now use a prior that is zero below p = 0.5 and a constant above p = 0.5. This corresponds to prior information that a majority of the Earth’s surface is water. Repeat each problem above and compare the inferences. Plot the posterior. What difference does the better prior make? If it helps, compare inferences (using both priors) to the true value p = 0.7.

#A

# Redetermining of Posterior distribution probability with the utilization of grid approximation approach and with a prior that is zero below p = 0.5 and a constant above p = 0.5.

Pstr_grid <- seq(from = 0, to = 1, length.out = 1000 )

Prior <- ifelse(Pstr_grid < 0.5, 0, 1)

# Therefore the new likelihood of the distribution under the assumption of 8 water in 15 tosses would be as under:

Likelihood <- dbinom(8, size = 15, prob = Pstr_grid)

Posterior <- Likelihood * Prior

Posterior <- Posterior / sum(Posterior)

# Plot construction for the posterior distribution.

Plot_Pstr_Distrbtn<- plot( x = Pstr_grid, y = Posterior, type = "l" )

Plot_Pstr_Distrbtn
## NULL
# B

## Determination and drawing 10,000 samples from the above calculated grid approximation:

Samples <- sample( Pstr_grid, prob = Posterior, size = 10000, replace = TRUE )

# Utilization of the sample for the calculation of the 90 % HPDI for P.

HPDI_Calcultn_90_Prcnt <- HPDI(Samples, prob = 0.90 )


HPDI_Calcultn_90_Prcnt
##      |0.9      0.9| 
## 0.5005005 0.7117117
# Hence the revised 90 % HPDI for p would be within [0.50 , 0.71], however our prior was broader and within [0.33 , 0.72] 


# C

Data <- rbinom(10000, size = 15, prob = Samples)

# Plotting the data histogram

Histgrm_Plot <- hist(Data)

Histgrm_Plot
## $breaks
##  [1]  2  3  4  5  6  7  8  9 10 11 12 13 14 15
## 
## $counts
##  [1]   53  136  344  653 1113 1593 1715 1730 1320  791  372  155   25
## 
## $density
##  [1] 0.0053 0.0136 0.0344 0.0653 0.1113 0.1593 0.1715 0.1730 0.1320 0.0791
## [11] 0.0372 0.0155 0.0025
## 
## $mids
##  [1]  2.5  3.5  4.5  5.5  6.5  7.5  8.5  9.5 10.5 11.5 12.5 13.5 14.5
## 
## $xname
## [1] "Data"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
# Determination for the probability of 8 water in 15 tosses.

Problty <- length(Data[Data == 8]) / length(Data)

Problty
## [1] 0.1593
# Hence the new probability of 8 water in 15 tosses would be 15.9 % and which is slightly greater as compared to our prior probability value of 14.5 %.


# D

# Determination of the probability of observing 6 water in 9 tosses.

Data_Rvsd_2 <- rbinom(10000, size = 9, prob = Samples)

# Plotting the new data histogram

Histgrm_Plot_Rvsd_2 <- hist(Data_Rvsd_2)

Histgrm_Plot_Rvsd_2
## $breaks
##  [1] 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0
## 
## $counts
##  [1]    3   47    0  254    0  794    0 1593    0 2300    0 2387    0 1688    0
## [16]  740    0  194
## 
## $density
##  [1] 0.0006 0.0094 0.0000 0.0508 0.0000 0.1588 0.0000 0.3186 0.0000 0.4600
## [11] 0.0000 0.4774 0.0000 0.3376 0.0000 0.1480 0.0000 0.0388
## 
## $mids
##  [1] 0.25 0.75 1.25 1.75 2.25 2.75 3.25 3.75 4.25 4.75 5.25 5.75 6.25 6.75 7.25
## [16] 7.75 8.25 8.75
## 
## $xname
## [1] "Data_Rvsd_2"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## Probability calculation of observing 6 water in 9 tosses

Problty_Rvsd <- length(Data_Rvsd_2[Data_Rvsd_2 == 6]) / length(Data_Rvsd_2)

Problty_Rvsd
## [1] 0.2387
## Therefore under our new condition the Probability calculation of observing 6 water in 9 tosses would be 0.24 or 24 percent and greater than our prior probability calculation of 17.2 percent.

3-6. Suppose you want to estimate the Earth’s proportion of water very precisely. Specifically, you want the 99% percentile interval of the posterior distribution of p to be only 0.05 wide. This means the distance between the upper and lower bound of the interval should be 0.05. How many times will you have to toss the globe to do this?

## Determination of precise Earth's water calculation by utilization of 99 % percentile interval of the Posterior distribution.


Pi_Width_Calculatn <- function(N, True_p) {
  Likelihood <- dbinom(round(N*True_p), size = N, prob = Pstr_grid)
  Posterior <- Likelihood * Prior
  Posterior <- Posterior / sum(Posterior)
  Samples <- sample( Pstr_grid, prob = Posterior, size = 10000, replace = TRUE )
  interval <- PI(Samples, prob = 0.99)
  names(interval) <- NULL
  diff(interval)
}

True_p <- 0.7

#1. PI width calculation when N = 100, yielded posterior distribution width of approx 0.23.

N <- 100

Pi_Width_Calculatn(N, True_p)
## [1] 0.2312362
#2. PI width calculation when N = 1000, yielded posterior distribution width of approx 0.075.

N <- 1000

Pi_Width_Calculatn(N, True_p)
## [1] 0.07507508
#3. PI width calculation when N = 2000, yielded posterior distribution width of approx 0.051.

N <- 2000

Pi_Width_Calculatn(N, True_p)
## [1] 0.05405405
#4. PI width calculation when N = 2200, yielded the expected posterior distribution width of approx 0.050.

N <- 2200

Pi_Width_Calculatn(N, True_p)
## [1] 0.05105105
#5. PI width calculation when N = 2300, yielded posterior distribution width of approx 0.049.

N <- 2300

Pi_Width_Calculatn(N, True_p)
## [1] 0.04904905
## Analysis Result Findings: On the test of data ranging from N = 10, to N = 2300, the toss times of N = 2200 was observed to be closest for the p to be only 0.05 wide, or the total times the globe must be tossed is 2300, for the the distance between the upper and lower bound of the interval to be around 0.05.