# BSPBSS

Bayesian Spatial Blind Source Separation via Thresholded Gaussian Process.

## Installation

Install the released version of BSPBSS from Github with:

``devtools::install_github("benwu233/BSPBSS")``

## A toy example

This is a basic example which shows you how to solve a common problem.

First we load the package and generate simulated images with a probabilistic ICA model:

``````library(BSPBSS)
set.seed(612)
sim = sim_2Dimage(length = 30, sigma = 5e-4, n = 30, smooth = 6)``````

The true source signals are three 2D geometric patterns (set `smooth=0` to generate patterns with sharp edges).

``levelplot2D(sim\$S,lim = c(-0.04,0.04), sim\$coords)``

which generate observed images such as

``levelplot2D(sim\$X[1:3,], lim = c(-0.12,0.12), sim\$coords)``

Then we generate initial values for mcmc,

``ini = init_bspbss(sim\$X, sim\$coords, q = 3, ker_par = c(0.1,50), num_eigen = 50)``

and run!

``````res = mcmc_bspbss(ini\$X,ini\$init,ini\$prior,ini\$kernel,n.iter=2000,n.burn_in=1000,thin=10,show_step=1000)
#> iter 1000 Sat Oct  1 22:12:33 2022
#>
#>  stepsize_zeta 0.0167947 accp_rate_zeta 0.34
#> iter 2000 Sat Oct  1 22:12:36 2022
#>
#>  stepsize_zeta 0.0167947 accp_rate_zeta 0.37``````

Then the results can be summarized by

``res_sum = sum_mcmc_bspbss(res, ini\$X, ini\$kernel, start = 101, end = 200, select_p = 0.5)``

and shown by

``levelplot2D(res_sum\$S, lim = c(-1.3,1.3), sim\$coords)``

For comparison, we show the estimated sources provided by informax ICA here.

``levelplot2D(ini\$init\$ICA_S, lim = c(-1.7,1.7), sim\$coords)``

We may overspecify the number of components and still obtain reasonable results.

``````ini = init_bspbss(sim\$X, sim\$coords, q = 5, ker_par = c(0.1,50), num_eigen = 50)
res = mcmc_bspbss(ini\$X,ini\$init,ini\$prior,ini\$kernel, n.iter=2000,n.burn_in=1000,thin=10,show_step=1000)
#> iter 1000 Sat Oct  1 22:12:40 2022
#>
#>  stepsize_zeta 0.0104649 accp_rate_zeta 0.35
#> iter 2000 Sat Oct  1 22:12:43 2022
#>
#>  stepsize_zeta 0.0104649 accp_rate_zeta 0.43
res_sum = sum_mcmc_bspbss(res, ini\$X, ini\$kernel, start = 101, end = 200, select_p = 0.5)
levelplot2D(res_sum\$S, lim = c(-1.3,1.3), sim\$coords)``````

## NIfTI data

Install a small pack contains simulated NIfTI images and MNI152 template.

``devtools::install_github("benwu233/SIMDATA")``

Load and take a glance at the data.

``````data_path = system.file("extdata",package="SIMDATA")
neurobase::ortho2(t1_3mm, y=sim_4d[,,,1])``````

Conduct BSPBSS.

``````X = pre_nii(sim_4d,mask)
ini = init_bspbss(X\$data, X\$coords, dens = 0.5, q = 3, ker_par = c(0.01, 120), num_eigen = 200)
res = mcmc_bspbss(ini\$X,ini\$init,ini\$prior,ini\$kernel,n.iter=3000,n.burn_in=2000,thin=10,show_step=1000,subsample_p=0.005)
#> iter 1000 Sat Oct  1 22:14:38 2022
#>
#>  stepsize_zeta 0.00296575 accp_rate_zeta 0.42
#> iter 2000 Sat Oct  1 22:16:28 2022
#>
#>  stepsize_zeta 0.00769239 accp_rate_zeta 0.42
#> iter 3000 Sat Oct  1 22:18:19 2022
#>
#>  stepsize_zeta 0.00769239 accp_rate_zeta 0.41
res_sum = sum_mcmc_bspbss(res, ini\$X, ini\$kernel, start = 201, end = 300, select_p = 0.5)``````

Output:

``ICs = output_nii(res_sum\$S,sim_4d,X\$coords,file=NULL,std=TRUE)``

Component 1:

``neurobase::ortho2(t1_3mm,y=ICs[,,,1])``

Component 2:

``neurobase::ortho2(t1_3mm,y=ICs[,,,2])``

Component 3:

``neurobase::ortho2(t1_3mm,y=ICs[,,,3])``