library(bliss)

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

1 One single functional covariate

1.1 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.

  set.seed(1)
  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

1.2 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).
  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)
#> Sample from the posterior distribution.
#> Gibbs Sampler: 
#>   Initialization.
#>   Determine the starting point.
#>   Start the Gibbs Sampler.
#>   10%
#>   20%
#>   30%
#>   40%
#>   50%
#>   60%
#>   70%
#>   80%
#>   90%
#>   100%
#>   Return the result.
#> Coefficient function: smooth estimate.
#> Coefficient function: Bliss estimate.
#> Compute the approximation of the posterior distribution.
#> Support estimation.
#> Compute the (log) densities of the posterior sample.
  
  # Structure of a Bliss object
  str(res_bliss)
#> List of 12
#>  $ alpha                 :List of 1
#>   ..$ : num [1:50] 0.03 0.036 0.0549 0.2557 0.6414 ...
#>  $ beta_posterior_density:List of 1
#>   ..$ :List of 4
#>   .. ..$ grid_t         : num [1:512] 0 0.00196 0.00391 0.00587 0.00783 ...
#>   .. ..$ grid_beta_t    : num [1:512] -3.73 -3.71 -3.7 -3.68 -3.67 ...
#>   .. ..$ density        : num [1:512, 1:512] 6.65e-18 1.00e-17 1.51e-17 2.27e-17 3.41e-17 ...
#>   .. ..$ new_beta_sample: num [1:800, 1:50] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ beta_sample           :List of 1
#>   ..$ : num [1:1001, 1:50] 0 0 1.39 1.41 0 ...
#>  $ Bliss_estimate        :List of 1
#>   ..$ : num [1:50, 1] 0 0 0 2.6 2.6 ...
#>  $ chains                : NULL
#>  $ chains_info           :List of 1
#>   ..$ :List of 2
#>   .. ..$ estimates   :List of 3
#>   .. .. ..$ mu_hat         : num 1
#>   .. .. ..$ sigma_sq_hat   : num 0.821
#>   .. .. ..$ Smooth_estimate: num [1:50] 0.0587 0.0456 0.0988 0.6937 1.9176 ...
#>   .. ..$ autocorr_lag: num [1:50, 1:3] -0.70929 -0.05692 0.11248 0.05391 0.00997 ...
#>   .. .. ..- attr(*, "dimnames")=List of 2
#>   .. .. .. ..$ : NULL
#>   .. .. .. ..$ : chr [1:3] "mu" "sigma_sq" "beta"
#>  $ data                  :List of 6
#>   ..$ Q     : num 1
#>   ..$ y     : num [1:100] 0.145 1.166 -0.224 2.462 0.118 ...
#>   ..$ x     :List of 1
#>   .. ..$ : num [1:100, 1:50] -1.733 -1.068 -1.283 -2.545 0.483 ...
#>   .. .. ..- attr(*, "scaled:center")= num [1:50] 1.148 1.066 0.96 0.887 0.788 ...
#>   ..$ betas :List of 1
#>   .. ..$ : num [1:50] 0 0 0 0 3 3 3 3 3 3 ...
#>   ..$ grids :List of 1
#>   .. ..$ : num [1:50] 0 0.0204 0.0408 0.0612 0.0816 ...
#>   ..$ x_save:List of 1
#>   .. ..$ : num [1:100, 1:50] -0.5847 0.0801 -0.1352 -1.3969 1.631 ...
#>  $ posterior_sample      :List of 3
#>   ..$ trace            : num [1:1001, 1:11] -0.0518 -1.3895 -0.2314 -0.3778 -0.373 ...
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : NULL
#>   .. .. ..$ : chr [1:11] "b_1" "b_2" "b_3" "m_1" ...
#>   ..$ param            :List of 6
#>   .. ..$ phi_l               :List of 1
#>   .. .. ..$ : num [1:49] 0.2263 0.1666 0.1227 0.0903 0.0665 ...
#>   .. ..$ K                   : num [1, 1] 3
#>   .. ..$ l_values_length     : num [1, 1] 49
#>   .. ..$ potential_intervals :List of 1
#>   .. .. ..$ : num [1:50, 1:49, 1:100] -1.537 -1.417 -1.224 -1.01 -0.859 ...
#>   .. ..$ grids               :List of 1
#>   .. .. ..$ : num [1:50] 0 0.0204 0.0408 0.0612 0.0816 ...
#>   .. ..$ normalization_values:List of 1
#>   .. .. ..$ : num [1:50, 1:49] 0.185 0.233 0.264 0.263 0.261 ...
#>   ..$ posterior_density: num [1:1001, 1:6] 0.00 1.53e-150 4.08e-49 5.06e-15 6.90e+01 ...
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : NULL
#>   .. .. ..$ : chr [1:6] "posterior density" "log posterior density" "likelihood" "log likelihood" ...
#>  $ Smooth_estimate       :List of 1
#>   ..$ : num [1:50, 1] 0.0587 0.0456 0.0988 0.6937 1.9176 ...
#>  $ support_estimate      :List of 1
#>   ..$ :'data.frame': 2 obs. of  2 variables:
#>   .. ..$ begin: num [1:2] 5 38
#>   .. ..$ end  : num [1:2] 30 50
#>  $ support_estimate_fct  :List of 1
#>   ..$ : num [1:50] 0 0 0 0 1 1 1 1 1 1 ...
#>  $ trace_sann            :List of 1
#>   ..$ : num [1:50001, 1:16] 0.601 0.601 0.601 0.601 -0.168 ...
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : NULL
#>   .. .. ..$ : chr [1:16] "b_1" "b_2" "b_3" "b_4" ...
#>  - attr(*, "class")= chr "bliss"

1.3 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.

  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")

  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)

2 Several functional covariates

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

2.1 Simulate the dataset

   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)

2.2 Apply the Bliss method

  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)

2.3 Plot the result

   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)

3 Session informations

#> R version 4.0.5 (2021-03-31)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Debian GNU/Linux 10 (buster)
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.3.5.so
#> 
#> locale:
#>  [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=C              
#>  [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
#>  [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] bliss_1.0.2
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.5        knitr_1.30        magrittr_2.0.1    kutils_1.70      
#>  [5] splines_4.0.5     MASS_7.3-53.1     statmod_1.4.35    xtable_1.8-4     
#>  [9] lattice_0.20-41   rlang_0.4.10      minqa_1.2.4       carData_3.0-4    
#> [13] stringr_1.4.0     plyr_1.8.6        tools_4.0.5       grid_4.0.5       
#> [17] nlme_3.1-152      xfun_0.19         htmltools_0.5.0   yaml_2.2.1       
#> [21] lme4_1.1-26       digest_0.6.27     Matrix_1.3-2      zip_2.1.1        
#> [25] nloptr_1.2.2.2    rockchalk_1.8.144 evaluate_0.14     rmarkdown_2.6    
#> [29] openxlsx_4.2.3    stringi_1.5.3     compiler_4.0.5    boot_1.3-27      
#> [33] foreign_0.8-81