BSPBSS-vignette

library(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)
#> Warning: Removed 177 rows containing missing values or values outside the scale range
#> (`geom_tile()`).

which generate observed images such as

levelplot2D(sim$X[1:3,], lim = c(-0.12,0.12), sim$coords)
#> Warning: Removed 177 rows containing missing values or values outside the scale range
#> (`geom_tile()`).

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=100)
#> iter 100 Sat Jul  4 23:34:05 2026
#> 
#> zeta0.122581 stepsize_zeta 0.00712258 accp_rate_zeta 0.37
#> iter 200 Sat Jul  4 23:34:05 2026
#> 
#> zeta0.225717 stepsize_zeta 0.00783484 accp_rate_zeta 0.39
#> iter 300 Sat Jul  4 23:34:06 2026
#> 
#> zeta0.212185 stepsize_zeta 0.00861832 accp_rate_zeta 0.38
#> iter 400 Sat Jul  4 23:34:06 2026
#> 
#> zeta0.169832 stepsize_zeta 0.00948015 accp_rate_zeta 0.48
#> iter 500 Sat Jul  4 23:34:06 2026
#> 
#> zeta0.185979 stepsize_zeta 0.0104282 accp_rate_zeta 0.38
#> iter 600 Sat Jul  4 23:34:06 2026
#> 
#> zeta0.163719 stepsize_zeta 0.011471 accp_rate_zeta 0.37
#> iter 700 Sat Jul  4 23:34:06 2026
#> 
#> zeta0.194074 stepsize_zeta 0.0126181 accp_rate_zeta 0.35
#> iter 800 Sat Jul  4 23:34:06 2026
#> 
#> zeta0.200016 stepsize_zeta 0.0138799 accp_rate_zeta 0.32
#> iter 900 Sat Jul  4 23:34:07 2026
#> 
#> zeta0.187286 stepsize_zeta 0.0152679 accp_rate_zeta 0.29
#> iter 1000 Sat Jul  4 23:34:07 2026
#> 
#> zeta0.184249 stepsize_zeta 0.0152679 accp_rate_zeta 0.21
#> iter 1100 Sat Jul  4 23:34:07 2026
#> 
#> zeta0.15656 stepsize_zeta 0.0152679 accp_rate_zeta 0.3
#> iter 1200 Sat Jul  4 23:34:07 2026
#> 
#> zeta0.183467 stepsize_zeta 0.0152679 accp_rate_zeta 0.32
#> iter 1300 Sat Jul  4 23:34:07 2026
#> 
#> zeta0.149271 stepsize_zeta 0.0152679 accp_rate_zeta 0.4
#> iter 1400 Sat Jul  4 23:34:08 2026
#> 
#> zeta0.161001 stepsize_zeta 0.0152679 accp_rate_zeta 0.27
#> iter 1500 Sat Jul  4 23:34:08 2026
#> 
#> zeta0.214896 stepsize_zeta 0.0152679 accp_rate_zeta 0.28
#> iter 1600 Sat Jul  4 23:34:08 2026
#> 
#> zeta0.216023 stepsize_zeta 0.0152679 accp_rate_zeta 0.24
#> iter 1700 Sat Jul  4 23:34:08 2026
#> 
#> zeta0.155854 stepsize_zeta 0.0152679 accp_rate_zeta 0.25
#> iter 1800 Sat Jul  4 23:34:08 2026
#> 
#> zeta0.186504 stepsize_zeta 0.0152679 accp_rate_zeta 0.27
#> iter 1900 Sat Jul  4 23:34:09 2026
#> 
#> zeta0.214563 stepsize_zeta 0.0152679 accp_rate_zeta 0.26
#> iter 2000 Sat Jul  4 23:34:09 2026
#> 
#> zeta0.19612 stepsize_zeta 0.0152679 accp_rate_zeta 0.24

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)
#> Warning: Removed 177 rows containing missing values or values outside the scale range
#> (`geom_tile()`).

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)
#> Warning: Removed 177 rows containing missing values or values outside the scale range
#> (`geom_tile()`).

- A toy example