# Examples demonstrate the VAR algorithm for variable selection in 
# heteroscedastic/homoscedastic linear regression.
# Reference: "Variational approximation for heteroscedastic models and matching 
#		  pursuit algorithms" by David Nott et al.
# Written by Minh-Ngoc Tran (ngoctm@nus.edu.sg)

library(MASS)
library(Matrix)
#========================================================#
# Examples for heteroscedastic regression
# the search for inclusion in the variance model is restricted
p = q = 10
tbeta =  c(5,-4,-3,2,0,rep(0,p-5))
talpha = c(0,-4, 3,0,0,rep(0,q-5))
n = 200; sig20 = 1;
Sig = matrix(rep(1,p*p),p)
for (ii in 1:p) for (jj in 1:p) if (ii!=jj) Sig[ii,jj] = 0.5^(abs(ii-jj))
XX = mvrnorm(n, rep(0,p), Sig)
XX = pnorm(XX)
sig2 = sig20*exp(XX%*%talpha)
sig = sqrt(sig2)
y = 2+XX%*%tbeta + sig*rnorm(n)
X = cbind(rep(1,n),XX)
fit = vsva(y,X)
true.beta = c(2,tbeta); true.alpha = c(log(sig20),talpha)
beta.hat = fit$mub; alpha.hat = fit$muv
cbind(true.beta,beta.hat,true.alpha,alpha.hat)

# the search for inclusion in the variance model isn't restricted
p = q = 10
tbeta =  c(5,-4,-3, 2,0,rep(0,p-5))
talpha = c(0, 0, 0,-4,3,rep(0,q-5))
n = 200; sig20 = 1;
Sig = matrix(rep(1,p*p),p)
for (ii in 1:p) for (jj in 1:p) if (ii!=jj) Sig[ii,jj] = 0.5^(abs(ii-jj))
XX = mvrnorm(n, rep(0,p), Sig)
XX = pnorm(XX)
sig2 = sig20*exp(XX%*%talpha)
sig = sqrt(sig2)
y = 2+XX%*%tbeta + sig*rnorm(n)
X = cbind(rep(1,n),XX)
fit = vsva(y,X,X,restrict=FALSE)
true.beta = c(2,tbeta); true.alpha = c(log(sig20),talpha)
beta.hat = fit$mub; alpha.hat = fit$muv
cbind(true.beta,beta.hat,true.alpha,alpha.hat)

# high-dimensional example
p = q = 500
tbeta =  c(5,-4,-3,2,0,rep(0,p-5))
talpha = c(0,-4, 3,0,0,rep(0,q-5))
n = 200; sig20 = 1;
Sig = matrix(rep(1,p*p),p)
for (ii in 1:p) for (jj in 1:p) if (ii!=jj) Sig[ii,jj] = 0.5^(abs(ii-jj))
XX = mvrnorm(n, rep(0,p), Sig)
XX = pnorm(XX)
sig2 = sig20*exp(XX%*%talpha)
sig = sqrt(sig2)
y = 2+XX%*%tbeta + sig*rnorm(n)
X = cbind(rep(1,n),XX)
fit = vsva(y,X)
true.beta = c(2,tbeta); true.alpha = c(log(sig20),talpha)
beta.hat = fit$mub; alpha.hat = fit$muv
cbind(true.beta,beta.hat,true.alpha,alpha.hat)


#========================================================#
# variable selection and estimation by the VAR algorithm for homoscedastic model
# small p
p = 10
n = 100; sigma = 1
tbeta = c(5,-4,3,-2,2,rep(0,p-5))
Sig = matrix(rep(1,p*p),p)
for (ii in 1:p) for (jj in 1:p) if (ii!=jj) Sig[ii,jj] = 0.5^(abs(ii-jj))
XX = mvrnorm(n, rep(0,p), Sig)
y = 2+XX%*%tbeta + sigma*rnorm(n)
X = cbind(rep(1,n),XX)
fit = vsva_homo(y,X)
true.beta = c(2,tbeta)
beta.hat = fit$mub
cbind(true.beta,beta.hat)

# large p
p = 1000
n = 100; sigma = 1
tbeta = c(5,-4,3,-2,2,rep(0,p-5))
Sig = diag(p)
XX = mvrnorm(n, rep(0,p), Sig)
y = 2+XX%*%tbeta + sigma*rnorm(n)
X = cbind(rep(1,n),XX)
fit = vsva_homo(y,X)
true.beta = c(2,tbeta)
beta.hat = fit$mub
cbind(true.beta,beta.hat)

