Introduction to BLiSS method
Paul-Marie Grollemund
2019-02-28
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 data1.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.0659 0.0759 0.1149 0.2977 0.6424 ...
#> $ 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.09 -3.08 -3.06 -3.04 -3.02 ...
#> .. ..$ density : num [1:512, 1:512] 7.99e-28 8.21e-28 8.43e-28 8.64e-28 8.83e-28 ...
#> .. ..$ 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 1.72 1.54 2.07 2.11 ...
#> $ Bliss_estimate :List of 1
#> ..$ : num [1:50, 1] 0 0 0 0 2.6 ...
#> $ chains : NULL
#> $ chains_info :List of 1
#> ..$ :List of 2
#> .. ..$ estimates :List of 3
#> .. .. ..$ mu_hat : num 0.836
#> .. .. ..$ sigma_sq_hat : num 0.173
#> .. .. ..$ Smooth_estimate: num [1:50] 0.142 0.164 0.268 0.802 1.913 ...
#> .. ..$ autocorr_lag: num [1:50, 1:3] -0.18229 0.02207 -0.01357 0.00183 -0.03369 ...
#> .. .. ..- attr(*, "dimnames")=List of 2
#> .. .. .. ..$ : NULL
#> .. .. .. ..$ : chr [1:3] "mu" "sigma_sq" "beta"
#> $ data :List of 6
#> ..$ Q : num 1
#> ..$ y : num [1:100] 1.422 -0.847 3.799 0.726 1.087 ...
#> ..$ x :List of 1
#> .. ..$ : num [1:100, 1:50] -0.854 0.093 3.415 0.355 -1.169 ...
#> .. .. ..- attr(*, "scaled:center")= num [1:50] 0.884 0.808 0.731 0.671 0.616 ...
#> ..$ 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.03 0.977 4.299 1.24 -0.284 ...
#> $ posterior_sample :List of 3
#> ..$ trace : num [1:1001, 1:11] -0.201 0.313 -0.113 -0.128 -0.112 ...
#> .. ..- 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.775 -0.781 -0.781 -0.774 -0.775 ...
#> .. ..$ 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.175 0.233 0.251 0.26 0.256 ...
#> ..$ posterior_density: num [1:1001, 1:6] 4.19e-85 1.31e-37 7.59e-10 6.80e-01 1.52e+09 ...
#> .. ..- 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.142 0.164 0.268 0.802 1.913 ...
#> $ support_estimate :List of 1
#> ..$ :'data.frame': 3 obs. of 2 variables:
#> .. ..$ begin: num [1:3] 5 25 40
#> .. ..$ end : num [1:3] 15 30 49
#> $ 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:13] 0.5675 0.5675 0.5675 -0.0539 -0.0539 ...
#> .. ..- attr(*, "dimnames")=List of 2
#> .. .. ..$ : NULL
#> .. .. ..$ : chr [1:13] "b_1" "b_2" "b_3" "m_1" ...
#> - 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.5.2 (2018-12-20)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 16.04.6 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/libblas/libblas.so.3.6.0
#> LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
#>
#> 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.0
#>
#> loaded via a namespace (and not attached):
#> [1] Rcpp_1.0.0 knitr_1.21 magrittr_1.5
#> [4] splines_3.5.2 kutils_1.60 MASS_7.3-51.1
#> [7] mnormt_1.5-5 lattice_0.20-38 pbivnorm_0.6.0
#> [10] xtable_1.8-3 minqa_1.2.4 carData_3.0-2
#> [13] stringr_1.3.1 plyr_1.8.4 tools_3.5.2
#> [16] grid_3.5.2 nlme_3.1-137 xfun_0.4
#> [19] htmltools_0.3.6 yaml_2.2.0 lme4_1.1-19
#> [22] digest_0.6.18 lavaan_0.6-3 Matrix_1.2-15
#> [25] zip_1.0.0 nloptr_1.2.1 rockchalk_1.8.129
#> [28] evaluate_0.12 rmarkdown_1.11 openxlsx_4.1.0
#> [31] stringi_1.2.4 compiler_3.5.2 stats4_3.5.2
#> [34] foreign_0.8-71