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

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

  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.148 0.182 0.288 0.591 0.845 ...
#>  $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] -6.17 -6.15 -6.13 -6.11 -6.09 ...
#>   .. ..$density : num [1:512, 1:512] 2.03e-22 3.29e-22 5.30e-22 8.54e-22 1.37e-21 ... #> .. ..$ 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 -3.41 1.88 3.27 ...
#>  $Bliss_estimate :List of 1 #> ..$ : num [1:50, 1] 0 0 0 2.57 2.57 ...
#>  $chains : NULL #>$ chains_info           :List of 1
#>   ..$:List of 2 #> .. ..$ estimates   :List of 3
#>   .. .. ..$mu_hat : num 1.04 #> .. .. ..$ sigma_sq_hat   : num 0.302
#>   .. .. ..$Smooth_estimate: num [1:50] 0.325 0.405 0.696 1.593 2.428 ... #> .. ..$ autocorr_lag: num [1:50, 1:3] -0.2915 -0.1568 -0.1449 -0.0898 -0.0187 ...
#>   .. .. ..- 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.924 1.837 0.91 1.684 -0.96 ... #> ..$ x     :List of 1
#>   .. ..$: num [1:100, 1:50] -0.204 -0.538 0.62 -0.298 -0.142 ... #> .. .. ..- attr(*, "scaled:center")= num [1:50] 0.905 0.821 0.738 0.668 0.6 ... #> ..$ 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.701 0.367 1.526 0.607 0.764 ... #>$ posterior_sample      :List of 3
#>   ..$trace : num [1:1001, 1:11] -0.271 1.414 1.655 0.503 1.024 ... #> .. ..- 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] -0.1813 -0.1336 -0.0463 0.106 0.2833 ...
#>   .. ..$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.165 0.217 0.233 0.241 0.235 ...
#>   ..$posterior_density: num [1:1001, 1:6] 0.00 7.91e-111 1.63e-62 2.18e-21 6.47e-27 ... #> .. ..- 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.325 0.405 0.696 1.593 2.428 ... #>$ support_estimate      :List of 1
#>   ..$:'data.frame': 3 obs. of 2 variables: #> .. ..$ begin: num [1:3] 4 25 42
#>   .. ..$end : num [1:3] 16 30 46 #>$ support_estimate_fct  :List of 1
#>   ..$: num [1:50] 0 0 0 1 1 1 1 1 1 1 ... #>$ trace_sann            :List of 1
#>   ..$: num [1:50001, 1:16] 0.532 0.532 0.532 0.532 0.532 ... #> .. ..- 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 3.6.1 (2019-07-05)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 7 x64 (build 7601) Service Pack 1
#>
#> Matrix products: default
#>
#> locale:
#> [1] LC_COLLATE=C                   LC_CTYPE=French_France.1252
#> [3] LC_MONETARY=French_France.1252 LC_NUMERIC=C
#> [5] LC_TIME=French_France.1252
#>
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base
#>
#> other attached packages:
#> [1] bliss_1.0.1
#>
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.4        knitr_1.25        magrittr_1.5
#>  [4] splines_3.6.1     kutils_1.69       MASS_7.3-51.4
#>  [7] mnormt_1.5-5      lattice_0.20-38   pbivnorm_0.6.0
#> [10] xtable_1.8-4      rlang_0.4.5       minqa_1.2.4
#> [13] carData_3.0-2     stringr_1.4.0     plyr_1.8.4
#> [16] tools_3.6.1       grid_3.6.1        nlme_3.1-141
#> [19] xfun_0.10         htmltools_0.4.0   yaml_2.2.1
#> [22] lme4_1.1-21       digest_0.6.25     lavaan_0.6-5
#> [25] Matrix_1.2-17     zip_2.0.4         nloptr_1.2.1
#> [28] rockchalk_1.8.144 evaluate_0.14     rmarkdown_1.16
#> [31] openxlsx_4.1.2    stringi_1.4.3     compiler_3.6.1
#> [34] boot_1.3-23       stats4_3.6.1      foreign_0.8-72