#####################################
# Bayesian Statistics: Laboratory 4 #
# 05/05/2023                        #
#####################################
# Some useful package 
library(rstudioapi)
setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) 

library(rstan)
library(ggplot2)
library(bayesplot)

##############################################
# Load the dataset and ecplore its structure #
##############################################
data <- readRDS('pest_data.RDS')
str(data)

# Select a subset of variables 
Svar <- c("complaints", "building_id", "traps", 
          "date", "live_in_super","age_of_building",
          "total_sq_foot", "average_tenant_age", 
          "monthly_average_rent", "floors", "month")
pest_data <- data[, Svar]

##############################################
# number of buildings 
N_buildings <- length(unique(pest_data$building_id))
N_buildings

# Some Exploratory analysis
ggplot(pest_data, aes(x = complaints)) + 
  geom_bar()

ggplot(pest_data, aes(x = traps, y = complaints)) + 
  geom_jitter()

# To do:
# At first, organise the data into a list 

# Then, open the file **simple_poisson_regression_void.stan** and complete:
  # the data block
  # the paramter block
  # the model block 

# Note: Prior specification - we will consider two weakly informative prior, 
# since we expect negative slope on traps and a positive intercept, 
# that is $$ \beta \sim \mathcal{N}(-0.25, 1) \quad \quad \alpha \sim \mathcal{N}(log(4), 1); $$

# Then save the file as **simple_poisson_regression.stan**

# Finally compile the model by means of the *stan_model*() function and 
# sample draws from the model by using the *sampling*() function

# Check the correctness by printing the posterior summary


## arrange data into a list; To complete
stan_dat_simple <- 
  
str(stan_dat_simple)

## compile the model (after wirting it )
comp_model_P <- stan_model('simple_poisson_regression.stan')

# Sampling (refresh = 0 avoid printing the trace)
fit_P_real_data <- sampling(comp_model_P, 
                            data = stan_dat_simple, 
                            refresh = 0)
# Posterior summary
print(fit_P_real_data, pars = c('alpha','beta'),
      probs = c(0.025, 0.5, 0.975))