---
title: "Introduction to BLiSS method"
author: "Paul-Marie Grollemund"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: yes
    number_sections: true
vignette: >
  %\VignetteIndexEntry{Introduction to BliSS method}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r,echo=TRUE,message=FALSE, warning=FALSE}
  library(bliss)
```

This vignette describes step by step how to use the BLiSS method, by:

* simulate data to illustrate the BLiSS model
* obtain a sample of the *a posteriori* distribution with a sampler of Gibbs
* plot the *a posteriori* distribution of the coefficient function and that of its support
* calculate the different Bayesian estimators and plot the results in a graph

# One single functional covariate

## Simulation of the data set

To obtain data, several characteristics must be specified to simulate the data, such as n (number of observations), p (number of measurement moments) curves $x_{i}(.)$), beta$\_$types (the shape of the coefficient function), and b$\_$inf and b$\_$sup (to define the domain of curves $x_{i}(.)$). Based on these parameters, we can use the *sim* function in the following code to simulate $x_{i}(.)$ curves, and real values $y_{i}$, from the functional linear regression model.


```{r ,eval=TRUE,include = TRUE}
  param <- list(                        # define the "param" to simulate data
                Q=1,                    # the number of functional covariate
                n=100,                  # n is the sample size and p is the
                p=c(50),                # number of time observations of the curves
                beta_types=c("smooth"), # define the shape of the "true" coefficient function
                grids_lim=list(c(0,1))) # Give the beginning and the end of the observation's domain of the functions.

  data <- sim(param) # Simulate the data
```

## Apply the Bliss method

To obtain an *a posteriori* sample, we use the Gibbs algorithm. We use the main function *fit$\_$Bliss* which calls sub-functions that allow us

* to sample the *a posteriori* distribution, 
* then to calculate the *a posteriori* distribution of the coefficient function, 
* to execute an optimization algorithm to calculate a constant estimate by pieces, 
* to calculate an estimation of the support and to calculate the density of the sample a posteriori (useful to calculate the BIC criterion). 

This main function requires a **param** list containing

* **iter**, the number of iterations of the Gibbs algorithm, 
* **burnin** heating time, 
* **K**, hyperparameter K of the Bliss model,
* **grids**, the moments of measurement of the curves $x_{i}(.)$, 
* **prior$\_$beta**, an argument specifying the distribution a priori of $\beta$ **Ridge$\_$Zellner** is only considered in this vignette
* and **phi$\_$l**, an argument specifying the *a priori* distribution of (only the **Gamma** option is considered in this vignette).

```{r ,eval=TRUE, include = TRUE}
  param <- list(            # define the required values of the Bliss method.
                iter=1e3,   # The number of iteration of the main numerical algorithm of Bliss.
                burnin=2e2, # The number of burnin iteration for the Gibbs Sampler
                K=c(3))     # The number of intervals of the beta

   
  res_bliss<-fit_Bliss(data=data,param=param,verbose=TRUE)
  
  # Structure of a Bliss object
  str(res_bliss)
```

## Plot the result

We give here the code to obtain representations of the *a posteriori* distribution. First, we give the code to obtain a posteriori probabilities $\alpha(t|D)$, relative to the support. Then, the *image$\_$Bliss* function is used to represent the subsequent distribution of the coefficient function.

```{r ,eval=TRUE, include = TRUE,fig.height=5,fig.width=7}
  param$ylim <- range(range(res_bliss$beta_posterior_density[[1]]$grid_beta_t),
                      c(-5,5))
  param$cols <- rev(heat.colors(100))
  image_Bliss(res_bliss$beta_posterior_density,param,q=1)
  lines(res_bliss$data$grids[[1]],res_bliss$Bliss_estimate[[1]],type="s",lwd=2)
  lines(res_bliss$data$grids[[1]],res_bliss$data$betas[[1]],col=2,lwd=2,type="s")
```


```{r ,eval=TRUE, include = TRUE,fig.height=5,fig.width=7}
  ylim <- range(range(res_bliss$Bliss_estimate[[1]]),
                range(res_bliss$Smooth_estimate[[1]]))
  plot_bliss(res_bliss$data$grids[[1]],
             res_bliss$Bliss_estimate[[1]],lwd=2)
  lines_bliss(res_bliss$data$grids[[1]],
             res_bliss$Smooth_estimate[[1]],lty=2)
```

# Several functional covariates

To avoid execution lengths, this section is not executed. Please, give it a try.

## Simulate the dataset

```{r ,eval=FALSE, include = TRUE}
   param <- list(Q=2,
                 n=300,
                 p=c(40,60),
                 beta_shapes=c("simple","smooth"),
                 grids_lim=list(c(0,1),c(0,2)))

  data <- sim(param)
```

## Apply the Bliss method
```{r ,eval=FALSE, include = TRUE}
  param <- list(       # define the required values of the Bliss method.
     iter=1e3,         # The number of iteration of the main numerical algorithm of Bliss.
     burnin=2e2,       # The number of burnin iteration for the Gibbs Sampler
     K=c(3,3))         # The number of intervals of the beta

  res_Bliss_mult <- fit_Bliss(data=data,param=param)
```

## Plot the result
```{r ,eval=FALSE, include = TRUE,fig.height=5,fig.width=7}
   q <- 1
   param$ylim <- range(range(res_Bliss_mult$beta_posterior_density[[q]]$grid_beta_t),
                       c(-5,5))
   param$cols <- rev(heat.colors(100))
   image_Bliss(res_Bliss_mult$beta_posterior_density,param,q=q)
   lines(res_Bliss_mult$data$grids[[q]],res_Bliss_mult$Bliss_estimate[[q]],type="s",lwd=2)
   lines(res_Bliss_mult$data$grids[[q]],res_Bliss_mult$data$betas[[q]],col=2,lwd=2,type="s")

  ylim <- range(range(res_Bliss_mult$Bliss_estimate[[q]]),
                 range(res_Bliss_mult$Smooth_estimate[[q]]))
   plot_bliss(res_Bliss_mult$data$grids[[q]],
              res_Bliss_mult$Bliss_estimate[[q]],lwd=2,ylim=ylim)
   lines(res_Bliss_mult$data$grids[[q]],
         res_Bliss_mult$Smooth_estimate[[q]],lty=2)


   q <- 2
   param$ylim <- range(range(res_Bliss_mult$beta_posterior_density[[q]]$grid_beta_t),
                       c(-5,5))
   param$cols <- rev(heat.colors(100))
   image_Bliss(res_Bliss_mult$beta_posterior_density,param,q=q)
   lines(res_Bliss_mult$data$grids[[q]],res_Bliss_mult$Bliss_estimate[[q]],type="s",lwd=2)
   lines(res_Bliss_mult$data$grids[[q]],res_Bliss_mult$data$betas[[q]],col=2,lwd=2,type="l")
   
   ylim <- range(range(res_Bliss_mult$Bliss_estimate[[q]]),
                 range(res_Bliss_mult$Smooth_estimate[[q]]))
   plot_bliss(res_Bliss_mult$data$grids[[q]],
              res_Bliss_mult$Bliss_estimate[[q]],lwd=2,ylim=ylim)
   lines(res_Bliss_mult$data$grids[[q]],
         res_Bliss_mult$Smooth_estimate[[q]],lty=2)
```

# Session informations
```{r session,echo=FALSE,message=FALSE, warning=FALSE}
  sessionInfo()
```
















