Introduction to BLiSS method"

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
  library(bliss)

This vignette describes step by step how to use the BLiSS method. Below, you can find the following implemented features:

One single functional covariate case

Simulate a data set

In order to simulate a proper dataset for Bliss application, some characteristics must be specified:

Based on these parameters, data can be simulated (curves $x_{i}(.)$ and real values $y_{i}$) from the functional linear regression model by using the sim function, as suggested in the following chunck.

  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

How to apply the Bliss method

In order to apply the Bliss method, the main function to use is fit$_$Bliss. This function provides the following outputs:

An important required argument of the previous function is param, which is a list containing:

Find below, an example of use of this function and a sketch of the structure of the returned object.

  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)

Graphical results

This section presents how to obtain main graphical results (posterior quantities) derived from the Bliss method.

Coefficient function

Considering Functional Linear Regression model (FLR), and the specific scalar-on-function case, the major model parameter to infer is the coefficient function $\beta(.)$. The following chunck shows how to plot the posterior 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)

Additionnaly to this plot, one could usually want to display a point estimate of the coefficient function (which is a function). By using the following code, you can access to: Bliss estimate, a piecewise constant version of the coefficient function, and the smooth estimate, the standard bayesian estimate of the coefficient function (standard means that it minimizes the posterior $L^2$-loss).

  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)

  # Bliss estimate
  lines(res_bliss$data$grids[[1]],res_bliss$Bliss_estimate[[1]],type="s",lwd=2) 

  # Smooth estimate
  lines(res_bliss$data$grids[[1]],res_bliss$Smooth_estimate[[1]],lty=2)

  # True coefficient function
  lines(res_bliss$data$grids[[1]],res_bliss$data$betas[[1]],col="purple",lwd=2,type="s")

The solid black line is Bliss estimate, the dashed black line is the smooth estimate and the solid purple line is the true coefficien function.

Below, the function lines_bliss is used to plot the step function without the vertical lines, linking differents steps.

  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)

  # Bliss estimate
  lines_bliss(res_bliss$data$grids[[1]],res_bliss$Bliss_estimate[[1]],lwd=3) 

  # Smooth estimate
  lines(res_bliss$data$grids[[1]],res_bliss$Smooth_estimate[[1]],lty=2)

  # True coefficient function
  lines_bliss(res_bliss$data$grids[[1]],res_bliss$data$betas[[1]],col="purple",lwd=3)

Support of coefficient function

According to the scientific problematic, one could aim to infer the coefficient function, but it is possible to alternatively focus only on the support of the coefficient function. In this case, the sign and the magniture of the coefficient function could be considered as nuisance parameters. Therefore, the Bliss method provides a specific estimation procedure for the support of the coefficient function (which relies on the posterior distribution of the coefficient function). It consists in deriving the posterior probabilities $\alpha(t|D)$, for each $t$ in the domain $\mathcal T$ of the functional data, which correspond to the probabilities (conditionnaly to the observed data) that the support of the coefficient function covers the time $t$.

To plot the posterior probabilities, you have to use the following code :

  plot(res_bliss$alpha[[1]],type="o",xlab="time",ylab="posterior probabilities")

From these posterior probabilities, the support estimate is derived by thresholding the probabilities. Without prior information guiding the estimation procedure, the default threshold is 0.5. The estimate support is then defined as the collection of time $t$ for which the posterior probability $\alpha(t|D) > 0.5$.

  plot(res_bliss$alpha[[1]],type="o",xlab="time",ylab="posterior probabilities")
  abline(h=0.5,col=2,lty=2)

  for(i in 1:nrow(res_bliss$support_estimate[[1]])){
  segments(res_bliss$support_estimate[[1]]$begin[i],0.05,
           res_bliss$support_estimate[[1]]$end[i],0.05,col="red"
           )
  points(res_bliss$support_estimate[[1]]$begin[i],0.05,col="red",pch="|",lwd=2)
  points(res_bliss$support_estimate[[1]]$end[i],0.05,col="red",pch="|",lwd=2)
  }

A resume of the support estimate is provided with:

res_bliss$support_estimate[[1]]

Multiple functional covariates

To avoid unnecesseray computational time, this section is not executed. You could figure out that the functions, objects and procedures are mostly similar to the previous one (single functional covariate case). The main differences are that:

Simulate a data set

   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)

How to 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)

Graphical results

   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

  sessionInfo()


Try the bliss package in your browser

Any scripts or data that you put into this service are public.

bliss documentation built on March 18, 2022, 5:46 p.m.