% testing the uBaL function
warning off;
addpath utilities;
%-------------------------------------------%
% test for variable selection and estimation in linear regression
%-------------------------------------------%
clear all
n = 200;
d = 8;
Sigma = ones(d,d);
for ii = 1:d
    for jj=1:d
        if ii~=jj 
            Sigma(ii,jj) = 0.5^(abs(ii-jj));
        end
    end
end
X = mvnrnd(zeros(1,d),Sigma,n);
tbeta = [3;1.5;0;0;2;0;0;0];
% Compare the following three ways of using the uBaL. It's working well,
% isn't it?
y = 5 + X*tbeta + randn(n,1);
fit = uBaL(X,y);
fit.lambda
fit.beta
y = 0 + X*tbeta + randn(n,1); % No intercept in the true model
fit = uBaL(X,y);
fit.lambda
fit.beta
fit = uBaL(center(X),center(y),[],[],[],0);
fit.lambda
fit.beta
%-------------------------------------------%
% test for group selection and estimation in linear regression
%-------------------------------------------%
clear all
n = 200;
ng = 4; % number of groups;
groups = cell(ng,1);
for i=1:ng
    groups{i}=[2*i-1,2*i];
end
tbeta = [2,2,0,0,4,4,0,0]';
d = length(tbeta);
Sig = ones(d,d);
for ii = 1:d
    for jj=1:d
        if ii~=jj 
            Sig(ii,jj) = 0.5^(abs(ii-jj));
        end
    end
end
X = mvnrnd(zeros(1,d),Sig,n);
y =  5 + X * tbeta + randn(n,1);
XX = center(X);
yy = center(y);
fit = uBaL(XX,yy,'normal',groups,[],0);
fit.lambda
fit.beta
%-------------------------------------------%
% test for variable selection and estimation in logistic regression
%-------------------------------------------%
clear all
n = 200;
d = 8;
rand('state',111344);
randn('state',155323);
Sigma = ones(d,d);
for ii = 1:d
    for jj=1:d
        if ii~=jj 
            Sigma(ii,jj) = 0.5^(abs(ii-jj));
        end
    end
end
X = mvnrnd(zeros(1,d),Sigma,n);
tbeta = [3;1.5;0;0;2;0;0;0];
mu=exp(X*tbeta)./(1+exp(X*tbeta)); %Note: no intercept in the true model
fre=10.*ones(n,1);
y=[binornd(fre,mu),fre];
fit = uBaL(X,y,'binomial',[],[],0);
fit.lambda
fit.beta
fit = uBaL(X,y,'binomial',[],[],1);
fit.lambda
fit.beta

%-------------------------------------------%
% test for group selection in logistic regression
%-------------------------------------------%
clear all
n = 200;
ng = 4; % number of groups;
groups = cell(ng,1);
for i=1:ng
    groups{i}=[2*i-1,2*i];
end
tbeta = [2,2,0,0,4,4,0,0]';
d = length(tbeta);
Sig = ones(d,d);
for ii = 1:d
    for jj=1:d
        if ii~=jj 
            Sig(ii,jj) = 0.5^(abs(ii-jj));
        end
    end
end
X = mvnrnd(zeros(1,d),Sig,n);
mu=exp(5+X*tbeta)./(1+exp(5+X*tbeta));
fre=10.*ones(n,1);
y=[binornd(fre,mu),fre];
fit = uBaL(X,y,'binomial',groups);
fit.lambda
fit.beta
%-------------------------------------------%
% test for factor selection in linear regression
%-------------------------------------------%
clear all;
n = 300;
nme = 5; % number of factors, each has three levels
Z = randn(n,nme); % latent variables
for i = 1:n
    for ime = 1:nme
        if Z(i,ime)<-.45 X(i,2*ime-1) = 1; end
        if Z(i,ime)>.45 X(i,2*ime) = 1; end
    end
end
for ime = 1:nme
    groups{ime} = [2*ime-1,2*ime];
end
MainEffect1 = X(:,groups{1});
MainEffect3 = X(:,groups{3});
XX = [MainEffect1,MainEffect3];
tbeta = [3,2,3,2]';
y = XX * tbeta + randn(n,1);
fit = uBaL(X,y,'normal',groups,[],0);
fit.lambda
fit.beta
%-------------------------------------------%
% test for main effect and interaction selection in linear regression
%-------------------------------------------%
clear all;
n = 300;
nme = 3; % number of main effects, each main effect has three levels
Z = randn(n,nme); % latent variables
for i = 1:n
    for ime = 1:nme
        if Z(i,ime)<-.45 X(i,2*ime-1) = 1; end
        if Z(i,ime)>.45 X(i,2*ime) = 1; end
    end
end
for ime = 1:nme
    groups{ime} = [2*ime-1,2*ime];
end
MainEffect1 = X(:,groups{1});
MainEffect2 = X(:,groups{2});
Interaction12 = [];
for i = 1:size(MainEffect1,2)
    for j = 1:size(MainEffect2,2)
        Interaction12 = [Interaction12,MainEffect1(:,i).*MainEffect2(:,j)];
    end
end
XX = [MainEffect1,MainEffect2,Interaction12];
tbeta = [3,2,3,2,1,1,1,1]';
y = XX * tbeta + randn(n,1);
fit = uBaL(X,y,'normal',groups,1,0);
fit.lambda
fit.beta
