## Conditional Probability Simulation
## Metropolis Algorithm
## Visiting the islands proportional to population size.

library(tidyverse)

# Distribution to explore, 6 islands with population proportional to value

# Theory

theta <- 1:6

theta_prob <- theta/sum(theta)

theta_df <- data_frame(theta = theta,
                       theta_prob = theta_prob)

ggplot(as.data.frame(theta_prob), aes(x = theta, y = theta_prob)) + geom_col()

# Using Metropolis to explore the distribution
# Note in the simulation we assume we do not know the distribution, just that the
# politician wants to visit all of the islands and spend a proportional amount of
# time on each island, proportional to the population size.

# theta is the poplation size, 1 to 6

B <- 100000
theta <- numeric(B)
current <- 1

theta[1] <- current

th.min <- 1; th.max <- 6

for ( i in 2:B){
  if (runif(1) < .5) {
    p.move <- (current - 1)/current
    ifelse( runif(1) < min( p.move, 1), 
                            current <- max( current - 1, th.min),
                            current)
    theta[i] <- current
  } else {
    p.move <- (current + 1)/current
    ifelse( runif(1) < min( p.move, 1),
                            current <- min( current + 1, th.max),
                            current)
    theta[i] <- current
  }
}

theta_df <- data.frame( time = 1:B,
                        theta = theta)

library(tidyverse)

# Plot the last 1000 sampled values

ggplot(tail(theta_df, 1000), aes(x=time, y = theta)) + geom_line() + coord_flip()

# The main idea with the Metropolis Algorithm is to investigate the distribution.
# Note that it has investigated the distribution very well.

ggplot(theta_df, aes(x = factor(theta))) + geom_bar() +  labs(x = "theta")

# Summary statistics computed on the sampled values from the distribution

mean(theta)
sd(theta)
median(theta)