friedman_demo

library(mvBayes)
library(BASS)

Generate Data

f<-function(x){
  10*sin(pi*x[,1]*x[,2])+20*(x[,3]-.5)^2+10*x[,4]+5*x[,5]
}

sigma<-1
nfunc = 50
tt = seq(0, 1, length.out = nfunc)  # functional variable grid
n = 500  # sample size
p = 9  # number of predictors other (only 4 are used)
X<-matrix(runif(n*p),n,p) # 9 non-functional variables, only first 4 matter
x<-cbind(rep(tt,each=n),kronecker(rep(1,nfunc),X)) # to get y
Y<-matrix(f(x),nrow=n)+rnorm(n*nfunc,0,sigma)

ntest = 1000
Xtest = matrix(runif(ntest * p), ntest, p)
x<-cbind(rep(tt,each=ntest),kronecker(rep(1,nfunc),Xtest)) # to get y
Ytest = matrix(f(x),nrow=ntest)+rnorm(ntest*nfunc,0,sigma)

Fit a multivariate BASS model

mod = mvBayes(
  bass,
  X,
  Y,
  nBasis=3
)
#> Starting mvBayes with 3 components, using 1 cores.
#> MCMC Start #-- May 28 08:05:11 --# nbasis: 0 
#> MCMC iteration 1000 #-- May 28 08:05:12 --# nbasis: 12 
#> MCMC iteration 2000 #-- May 28 08:05:12 --# nbasis: 13 
#> MCMC iteration 3000 #-- May 28 08:05:12 --# nbasis: 13 
#> MCMC iteration 4000 #-- May 28 08:05:12 --# nbasis: 13 
#> MCMC iteration 5000 #-- May 28 08:05:12 --# nbasis: 11 
#> MCMC iteration 6000 #-- May 28 08:05:12 --# nbasis: 11 
#> MCMC iteration 7000 #-- May 28 08:05:12 --# nbasis: 11 
#> MCMC iteration 8000 #-- May 28 08:05:12 --# nbasis: 11 
#> MCMC iteration 9000 #-- May 28 08:05:12 --# nbasis: 11 
#> MCMC iteration 10000 #-- May 28 08:05:13 --# nbasis: 11 
#> MCMC Start #-- May 28 08:05:13 --# nbasis: 0 
#> MCMC iteration 1000 #-- May 28 08:05:13 --# nbasis: 8 
#> MCMC iteration 2000 #-- May 28 08:05:13 --# nbasis: 8 
#> MCMC iteration 3000 #-- May 28 08:05:13 --# nbasis: 9 
#> MCMC iteration 4000 #-- May 28 08:05:13 --# nbasis: 9 
#> MCMC iteration 5000 #-- May 28 08:05:13 --# nbasis: 9 
#> MCMC iteration 6000 #-- May 28 08:05:13 --# nbasis: 8 
#> MCMC iteration 7000 #-- May 28 08:05:13 --# nbasis: 8 
#> MCMC iteration 8000 #-- May 28 08:05:13 --# nbasis: 8 
#> MCMC iteration 9000 #-- May 28 08:05:14 --# nbasis: 8 
#> MCMC iteration 10000 #-- May 28 08:05:14 --# nbasis: 8 
#> MCMC Start #-- May 28 08:05:14 --# nbasis: 0 
#> MCMC iteration 1000 #-- May 28 08:05:14 --# nbasis: 8 
#> MCMC iteration 2000 #-- May 28 08:05:14 --# nbasis: 7 
#> MCMC iteration 3000 #-- May 28 08:05:14 --# nbasis: 7 
#> MCMC iteration 4000 #-- May 28 08:05:14 --# nbasis: 8 
#> MCMC iteration 5000 #-- May 28 08:05:14 --# nbasis: 7 
#> MCMC iteration 6000 #-- May 28 08:05:14 --# nbasis: 7 
#> MCMC iteration 7000 #-- May 28 08:05:14 --# nbasis: 7 
#> MCMC iteration 8000 #-- May 28 08:05:14 --# nbasis: 7 
#> MCMC iteration 9000 #-- May 28 08:05:15 --# nbasis: 7 
#> MCMC iteration 10000 #-- May 28 08:05:15 --# nbasis: 7
#> Generating 'samples' attribute, since it was absent in 'bmList[[1]]'
#> Approximating 'residSD', since 'residSDExtract' is None
plot(mod)
#> Warning in predict.mvBayes(x, Xtest, returnPostCoefs = TRUE, idxSamples =
#> idxSamples): 'idxSamples' is not an argument of the bayesModel predict
#> function...setting idxSamples='default'

Plot PCA Decomposition

plot(mod$basisInfo, idxMV = tt, xlabel = "tt")

Plot Traceplots

traceplot(mod)

Evaluate Training Fit

plot(mod, idxMV = tt, xlabel = "tt")
#> Warning in predict.mvBayes(x, Xtest, returnPostCoefs = TRUE, idxSamples =
#> idxSamples): 'idxSamples' is not an argument of the bayesModel predict
#> function...setting idxSamples='default'

Evaluate Test Fit

plot(mod, Xtest = Xtest, Ytest = Ytest, idxMV = tt, xlabel = "tt")
#> Warning in predict.mvBayes(x, Xtest, returnPostCoefs = TRUE, idxSamples =
#> idxSamples): 'idxSamples' is not an argument of the bayesModel predict
#> function...setting idxSamples='default'

Compute Sensitivity

modSensitivity = mvSobol(mod)
#> Start ##------ Thu May 28 08:05:20 2026 ------## 
#> Integrating ##------ Thu May 28 08:05:20 2026 ------## 
#> Shuffling ##------ Thu May 28 08:05:20 2026 ------## 
#> Finish ##------ Thu May 28 08:05:20 2026 ------##
#> Warning in mvSobol(mod): 'BASS' has not implemented totalOrder
plot(modSensitivity, idxMV = tt, xlabel = "tt")
#> Warning in plot.sobol(modSensitivity, idxMV = tt, xlabel = "tt"): 'x' does not
#> have totalOrderSobol

Compute Posterior Samples and Means

# All posterior predictive samples
Ytest_postSamples = predict(mod, Xtest)
# Posterior predictive mean
Ytest_postMean = apply(Ytest_postSamples, 2, mean)
# single posterior predictive sample (from MCMC iteration #429)
Ytest_postSample429 = predict(mod, Xtest, idxSamples = 429)
#> Warning in predict.mvBayes(mod, Xtest, idxSamples = 429): 'idxSamples' is not
#> an argument of the bayesModel predict function...setting idxSamples='default'