Gibbs_alg.Rmd

---
title: "Gibbs_alg"
author: "Renjie Wei"
date: '2022-05-04'
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
```{r}
library(MASS)
library(tidyverse)
library(extraDistr)
library(matrixsampling)
library(parallel)
library(foreach)
library(doParallel)
```

```{r}
load("./dt_long.RData")
dt_wind =
  dt_long %>% 
    mutate(intercept = 1) %>% 
    dplyr::select(c("ID", "Wind.kt", "intercept", "Wind_prev", "Lat_change", "Long_change", "Wind_change")) %>%
    drop_na()

# 699 after drop_na
IDs <- unique(dt_wind$ID)

# 684 observations >=5
ID_in <- NULL
for (i in 1:length(IDs)) {
    dt_temp =  dt_wind[which(dt_wind$ID == IDs[i]),][,2:7]
    m = nrow(dt_temp)
    # need to drop observations less than 5
    if (m <= 3){
        next
    }else{
        ID_in = c(ID_in,  IDs[i])
    }
}

dt_mtx_full = dt_wind %>% filter(ID %in% ID_in) %>% nest() %>% pull(data)

dt_mtx_train <- list(NULL)
dt_mtx_test <- list(NULL)

# partitioning dataset based on id 
# ##########dont do it!###########
#set.seed(20228160)
#ID_train = sample(ID_in, ceiling(0.8 * length(ID_in)))
#ID_test = NULL
#
#for(i in 1:length(ID_in)){
#    if(ID_in[i] %in% ID_train){
#       next
#    }else{
#        ID_test = c(ID_test, ID_in[i]) 
#    }
#}
#
#dt_mtx_train = dt_wind %>% filter(ID %in% ID_train) %>% nest() %>% pull(data)
#dt_mtx_test = dt_wind %>% filter(ID %in% ID_test) %>% nest() %>% pull(data)

#for (i in 1:length(ID_train)) {
#    dt_temp = dt_wind[which(dt_wind$ID == ID_train[i]),][,2:7]
#    dt_temp = dt_temp %>% as.matrix()
#    dt_mtx_train[[i]] = dt_temp
#}
#
#for (i in 1:length(ID_test)) {
#    dt_temp = dt_wind[which(dt_wind$ID == ID_test[i]),][,2:7]
#    dt_temp = dt_temp %>% as.matrix()
#    dt_mtx_test[[i]] = dt_temp
#}

# right way!
# another way to partitioning dataset, each i set 80% of observations to training set
dt_train = list(NULL)
dt_test = list(NULL)
set.seed(20228160)
for(i in 1:length(dt_mtx_full)){
    dt_temp = dt_mtx_full[[i]]
    rowTrain = sample(nrow(dt_temp), ceiling(0.8 * nrow(dt_temp)))
    train_ = dt_temp[rowTrain,]
    test_ = dt_temp[-rowTrain,]
    dt_train[[i]] = train_
    dt_test[[i]] = test_
}
```

```{r posterior_samp}
# posterior of B
# each beta_i 1*5
# B n*5

post.B <- function(dt, muvec, sigma2, Sigma){
   n = length(dt)
   B = NULL
   RSS = NULL
   m = NULL
   for (i in 1:n){
       # stuffs to define the distribution
       X = as.matrix(dt[[i]][,-1]) 
       y = as.vector(dt[[i]][,1])
       V = sigma2^(-1) * t(X) %*% X + solve(Sigma)
       M = sigma2^(-1) * t(y) %*% X + muvec %*% solve(Sigma)
       mean_bi = solve(V) %*% t(M) # we need a 5*1
       vcov_bi = solve(V)
       bi = mvrnorm(1, mu = mean_bi, Sigma = vcov_bi)
       B = rbind(B, bi)
       # calculate RSS,m_i for sigma2
       RSS = rbind(RSS, sum((y - X %*% bi)^2))
       m = rbind(m, nrow(X))
   }
   return(list(B = B, RSS = RSS, m = m))
}

# test passed
#testB = post.B(dt = dt_mtx, muvec = rep(0,5), sigma2 = 2, Sigma = diag(2,5,5))


# sampling sigma2 from a inverse gamma, need sum(m_i) and SSR
post.sigma2 <- function(m, RSS){
    alpha_ = sum(m)/2
    beta_ = sum(RSS)/2
    sigma2 = rinvgamma(1, alpha = alpha_, beta = beta_)
    return(sigma2)
}

# test passed
#testsigma2 = post.sigma2(testB$m, testB$RSS)

# sampling Sigma from inverse wishart distribution
post.Sigma <- function(B, muvec){
    n = nrow(B)
    S.matrix = diag(1, 5, 5)
    for (i in 1:n){
        beta_i = B[i,]
        S.matrix = S.matrix + (beta_i-muvec) %*% t(beta_i-muvec)
    }
    v = n + 5 + 1
    Sigma = rinvwishart(1, nu = v, Omega = S.matrix, checkSymmetry = F)
    return(Sigma[,,1])
}

# test passed
#testSigma = post.Sigma(testB$B, rep(1,5))

# sampling muvec from multivariate normal distribution
post.mu <- function(B, Sigma){
    n = nrow(B)
    mean_mu = colMeans(B)
    vcov_mu = 1/n * Sigma
    muvec = mvrnorm(1, mean_mu, vcov_mu)
    return(muvec)
}

# test passed
#testmuvec = post.mu(testB$B, testSigma)
```

```{r MCMC_Gibbs}
MCMC.Gibbs <- function(dt, init.B, init.muvec, init.sigma2, init.Sigma, max.iter = 1e4){
    B.res = list()
    sigma2.res = list()
    Sigma.res = list()
    muvec.res = list()
    
    # initialize B
    init.Bstuff = post.B(dt, init.muvec, init.sigma2, init.Sigma)
    
    cur.stuff = NULL
    cur = list(B = init.B, sigma2 = init.sigma2, muvec = init.muvec, Sigma = init.Sigma ,RSS = init.Bstuff$RSS)
    
    B.res[[1]] = init.B
    sigma2.res[[1]] = cur$sigma2
    Sigma.res[[1]] = cur$Sigma
    muvec.res[[1]] = cur$muvec
    
    prev.stuff = cur.stuff
    prev = cur
    
    # start iteration
    for (i in 2:max.iter) {
        cur.stuff = post.B(dt, prev$muvec, prev$sigma2, prev$Sigma)
        sigma2 = post.sigma2(cur.stuff$m, cur.stuff$RSS)
        Sigma = post.Sigma(cur.stuff$B, prev$muvec)
        muvec = post.mu(cur.stuff$B, Sigma)
        
        # update current parameters
        cur = list(B = cur.stuff$B, sigma2 = sigma2, muvec = muvec, Sigma = Sigma, RSS = cur.stuff$RSS)
        prev.stuff = cur.stuff
        prev = cur
        
        # append to final result
        B.res[[i]] = cur$B
        sigma2.res[[i]] = cur$sigma2
        Sigma.res[[i]] = cur$Sigma
        muvec.res[[i]] = cur$muvec
        
    }
    return(list(
        B = B.res,
        sigma2 = sigma2.res,
        Sigma = Sigma.res,
        muvec = muvec.res
    ))
}
```

```{r ols_beta}
ols.beta.res = NULL
ols.residual.res = NULL
for(i in 1:length(dt_train)){
    lm.mdl = lm(Wind.kt~ Wind_prev+ Lat_change+ Long_change+ Wind_change, data = dt_train[[i]])
    ols.beta.res = rbind(ols.beta.res, coef(lm.mdl))
    ols.residual = mean((lm.mdl$residuals)^2)
    ols.residual.res = rbind(ols.residual.res, ols.residual)
}
ols.B = as.matrix(ols.beta.res)
```



```{r run_MCMC, cache=T, eval=FALSE}
# train set based on ID
#MCMC.res = MCMC.Gibbs(dt_mtx_train, init.muvec = c(1,2,3,4,5), init.sigma2 = 5, init.Sigma = #diag(c(1,2,3,4,5)))
#save(MCMC.res,file = "MCMC.res.RData")
#finished

# train set within each hurricane
mu0 = colMeans(ols.B,na.rm = T)
sigmasq0 = mean(ols.residual.res)

# run 10 cycles
#test.MCMC.res = MCMC.Gibbs(dt_train, init.B = ols.B  ,init.muvec = mu0, init.sigma2 = sigmasq0, init.Sigma = diag(c(1,2,3,4,5)), max.iter = 10)

MCMC.res.2 = MCMC.Gibbs(dt_train, init.B = ols.B, init.muvec = mu0, init.sigma2 = sigmasq0, init.Sigma = diag(c(1,2,3,4,5)))
save(MCMC.res.2,file = "MCMC.res.2.RData")

bad.B = matrix(seq(1:(5*697)), 697)
bad.mu = colMeans(bad.B)
MCMC.res.bad =  MCMC.Gibbs(dt_train, init.B = bad.B, init.muvec = bad.mu, init.sigma2 = sigmasq0*1000, init.Sigma = diag(c(1000,50,1000,500,100)), max.iter = 10000)
save(MCMC.res.bad,file = "MCMC.res.bad.RData")
```

```{r}
#save(dt_mtx_train,file = "dt_mtx_train.RData")
#save(dt_mtx_test,file = "dt_mtx_test.RData")
save(dt_train,file = "dt_train.RData")
save(dt_test,file = "dt_test.RData")
```

```{r}
save(ID_in, file = "ID_in.RData")
```

```{r try_bad_starting_value}
# do we need to do that?

```


```{r ols_rss_train}
ols_rss_train = NULL
for(i in 1:length(dt_train)){
    lm.mdl = lm(Wind.kt~ Wind_prev+ Lat_change+ Long_change+ Wind_change, data = dt_train[[i]])
    ols_rss_train = rbind(ols_rss_train , sum((predict(lm.mdl, newdata = dt_train[[i]]) - dt_train[[i]][,"Wind.kt"])^2))
}
```

```{r mcmc_rss_train}
mcmc_rss_train = NULL
for (i in 1:length(dt_train)){
    betas = as.matrix(beta.res.postmean[i,1:5])
    mcmc_rss_train = rbind(mcmc_rss_train, sum((as.matrix(dt_train[[i]][2:6]) %*% t(betas) - dt_train[[i]][,1])^2))
}
```

```{r}
colMeans(mcmc_rss_train)
colMeans(ols_rss_train)
```