function retSD = SimuSGLMData(inFamily,inLink,modid,n,p)
%function retSD = SimuSGLMData(inFamily,inLink,modid,n,p)
%
% Simulate data from a sparse generalized linear model (GLM).
%
% Model: g(E[Y]) = mu + X*beta
%
% Input:
%   inFamily: family of the distribution;
%   inLink: link function of the generalized linear model;
%   modid: index of the models;
%   n: number of observation in each dataset;
%   p: number of variables in X.
%
% Output:
%   retSD = [Y,X,mu,beta]
%
% Dabao Zhang     November 7, 2011
%
% Revised by:
%       Dabao Zhang, May 21, 2013 -- Family & Link.
%
% Copyright (c) 2011-12 by Dabao Zhang.

error(nargchk(3, 5, nargin));

if( isempty(inLink) )
    switch inFamily
        case 'normal'
            inLink = 'identity'; 
        case 'binomial'
            inLink = 'logit';
        case 'poisson'
            inLink = 'log';
        case 'gamma'
            inLink = -1;
        case 'inverse gaussian'
            inLink = -2;
        case 'quasi' % Incomplete
            inLink = 'identity';
        otherwise
            error('  SimuSGLMData: The family is not supported...');
    end
end

switch inFamily
    case 'binomial'
        retSD = SimuSBinomialData(inLink,modid,n,p);
    case 'normal'
        retSD = SimuSGaussianData(inLink,modid,n,p);
    case 'poisson'
        retSD = SimuSPoissonData(inLink,modid,n,p);
    case 'gamma'
        retSD = SimuSGammaData(inLink,modid,n,p);
    case 'inverse gaussian'
        retSD = SimuSInvGaussianData(inLink,modid,n,p);
    case 'quasi'
        retSD = SimuSQuasiData(inLink,modid,n,p);
    otherwise
        error('  SimuSGLMData: The family is not supported...');
end



%-------------------------------------------------------------------------
function retSD = SimuSBinomialData(inLink,modid,n,p)
%function retSD = SimuSBinomialData(inLink,modid,n,p)
%
% Simulate data from a sparse binomial model.
%
% Model: logit(E[Y]) = mu + X*beta
%
% Input:
%   inLink: link function of the generalized linear model;
%   modid: index of the model;
%   n: number of observation in each dataset;
%   p: number of variables in X.
%
% Output:
%   retSD: = [Y,X,mu,beta,p]
%
% Dabao Zhang     November 7, 2011
%
% Revised by:
%
% Copyright (c) 2011-12 by Dabao Zhang.

retSD.X = zeros(n,p);
retSD.Y = zeros(n,1);
retSD.mu = 0;
retSD.beta = zeros(p,1);

switch modid
    case {1,2}
        retSD.beta(1:10) = 2;  retSD.beta(101:110) = 1;
        if( modid==1 ) rho = 0.9; else rho = 0.5; end
        sigma = sqrt(1-rho^2);

        nBlk = max(p/100,1);
        tmpBLen = min(p,100);
        % Here it is assumed that p = nBlk*100 or p<100
        for j=1:nBlk
            idS = 100*(j-1)+1;
            retSD.X(:,idS) = normrnd(0,1,n,1);
            for i=2:tmpBLen
                retSD.X(:,idS+i-1) = rho*retSD.X(:,idS+i-2)...
                                     +sigma*normrnd(0,1,n,1);
            end
        end
        eta = retSD.mu+retSD.X*retSD.beta;
        
    case 3
        retSD.beta(1:30) = 1.5;
        Z = normrnd(0,1,n,3);
        retSD.X(:,1:10) = Z(:,1)*ones(1,10)+0.1*normrnd(0,1,n,10);
        retSD.X(:,11:20) = Z(:,2)*ones(1,10)+0.1*normrnd(0,1,n,10);
        retSD.X(:,21:30) = Z(:,3)*ones(1,10)+0.1*normrnd(0,1,n,10);
        retSD.X(:,31:p) = normrnd(0,1,n,p-30);
        eta = retSD.mu+retSD.X*retSD.beta;
        
    case 4
        retSD.beta(1:5) = 1/10;   retSD.beta(6:10) = -1/10;
        retSD.beta(11:15) = 1/5; retSD.beta(16:20) = -1/5;
        retSD.beta(21:30) = 1/10;
        Z = normrnd(0,1,n,3);
        eta = retSD.mu+Z(:,1)+2*Z(:,2)+Z(:,3);
        retSD.X(:,1:5) = Z(:,1)*ones(1,5)+normrnd(0,1,n,5);
        retSD.X(:,6:10) = -Z(:,1)*ones(1,5)+normrnd(0,1,n,5);
        retSD.X(:,11:15) = Z(:,2)*ones(1,5)+normrnd(0,1,n,5);
        retSD.X(:,16:20) = -Z(:,2)*ones(1,5)+normrnd(0,1,n,5);
        retSD.X(:,21:30) = Z(:,3)*ones(1,10)+normrnd(0,1,n,10);
        retSD.X(:,31:p) = normrnd(0,1,n,p-30);

    otherwise
        error('This bionomial model is NOT implemented!');
end
retSD.Y = binornd(1,ginv(eta,inLink));



%-------------------------------------------------------------------------
function retSD = SimuSGaussianData(inLink,modid,n,p)
%function retSD = SimuSGaussianData(inLink,modid,n,p)
%
% Simulate data from a sparse Gaussian model.
%
% Model: E[Y] = mu + X*beta
%
% Input:
%   inLink: link function of the generalized linear model;
%   modid: index of the model;
%   n: number of observation in each dataset;
%   p: number of variables in X.
%
% Output:
%   retSD: = [Y,X,mu,beta,p]
%
% Dabao Zhang     November 7, 2011
%
% Revised by:
%
% Copyright (c) 2011-12 by Dabao Zhang.

retSD.X = zeros(n,p);
if( modid~=5 )
    retSD.mu = 0;
    retSD.beta = zeros(p,1);
    retSD.Y = zeros(n,1);
else
    retSD.mu = zeros(1,5);
    retSD.beta = zeros(p,5);
    retSD.Y = zeros(n,5);
end

% Case 1,2,3
switch modid
    case 1
        rho = 0.9;
    
    case 2
        rho = 0.5;

    case 5
        rho = 0.3;
end

switch modid
    case {1,2} % High/Mild Correlation
        retSD.beta(1:110) = [2*ones(10,1); zeros(90,1); ones(10,1)];
        sigma = sqrt(1-rho^2);

        nBlk = p/100;
        tmpBLen = min(p,100);
        % Here it is assumed that p = nBlk*100
        for j=1:nBlk
            idS = 100*(j-1)+1;
            retSD.X(:,idS) = normrnd(0,1,n,1);
            for k=2:tmpBLen
                retSD.X(:,idS+k-1) = rho*retSD.X(:,idS+k-2)...
                                     +sigma*normrnd(0,1,n,1);
            end
        end
        retSD.Y = retSD.mu+retSD.X*beta+normrnd(0,1,n,1);
        
    case 3 % Clusering
        retSD.beta(1:30) = 1.5;
        sigma=15;

        Z = normrnd(0,1,n,3);
        retSD.X(:,1:10) = Z(:,1)*ones(1,10)+0.1*normrnd(0,1,n,10);
        retSD.X(:,11:20) = Z(:,2)*ones(1,10)+0.1*normrnd(0,1,n,10);
        retSD.X(:,21:30) = Z(:,3)*ones(1,10)+0.1*normrnd(0,1,n,10);
        retSD.X(:,31:p) = normrnd(0,1,n,p-30);
        retSD.Y = retSD.mu+retSD.X*beta+sigma*normrnd(0,1,n,1);
        
    case 4 % Measurement Error Model
        %retSD.beta(1:3) = [1; 2; 1];
        retSD.beta(1:5) = 1/10;   retSD.beta(6:10) = -1/10;
        retSD.beta(11:15) = 1/5; retSD.beta(16:20) = -1/5;
        retSD.beta(21:30) = 1/10;

        Z = normrnd(0,1,n,3);
        retSD.Y = mu+Z(:,1)+2*Z(:,2)+Z(:,3)+normrnd(0,1,n,1);
        retSD.X(:,1:5) = Z(:,1)*ones(1,5)+normrnd(0,1,n,5);
        retSD.X(:,6:10) = -Z(:,1)*ones(1,5)+normrnd(0,1,n,5);
        retSD.X(:,11:15) = Z(:,2)*ones(1,5)+normrnd(0,1,n,5);
        retSD.X(:,16:20) = -Z(:,2)*ones(1,5)+normrnd(0,1,n,5);
        retSD.X(:,21:30) = Z(:,3)*ones(1,10)+normrnd(0,1,n,10);
        retSD.X(:,31:p) = normrnd(0,1,n,p-30);
        
    case 5 % Latent-Variable Model
        retSD.beta([50 150 250 350 450 550],1) = 2;
        retSD.beta([51 153 256 359 467 583],1) = -2;
        retSD.beta([50 51 150 153 250 256 350 359 450 467 550 583],2) = 2;
        retSD.beta([50 51 150 153 250 256 350 359 450 467 550 583],3) = -2;
        retSD.beta([50 150 250 350 450 550],[4 5]) = 3;
        retSD.beta([51 153 256 359 467 583],4) = 1.5;
        retSD.beta([51 153 256 359 467 583],5) = -1.5;
        sigma = sqrt(1-rho^2);

        nBlk = p/100;
        % Here it is assumed that p = nBlk*100
        for j=1:nBlk
            idS = 100*(j-1)+1;
            retSD.X(:,idS) = normrnd(0,1,n,1);
            for i=2:100
                retSD.X(:,idS+i-1) = rho*retSD.X(:,idS+i-2)...
                                     +sigma*normrnd(0,1,n,1);
            end
        end
        retSD.Y = ones(n,1)*retSD.mu+retSD.X*beta+normrnd(0,1,n,5);
        
    otherwise
        error('This model is NOT supported!');
    
end



%-------------------------------------------------------------------------
function retSD = SimuSPoissonData(inLink,modid,n,p)
%function retSD = SimuPoissonData(inLink,modid,n,p)
%
% Simulate data from a sparse Poisson model.
%
% Model: log(E[Y]) = mu + X*beta
%
% Input:
%   inLink: link function of the generalized linear model;
%   modid: index of the model;
%   n: number of observation in each dataset;
%   p: number of variables in X.
%
% Output:
%   retSD: = [Y,X,mu,beta,p]
%
% Dabao Zhang     November 4, 2011
%
% Revised by:
%
% Copyright (c) 2011-12 by Dabao Zhang.

error('This Poisson model is NOT implemented!');

retSD.X = zeros(n,p);
retSD.Y = zeros(n,1);
retSD.mu = 0;
retSD.beta = zeros(p,1);


          
%-------------------------------------------------------------------------
function retSD = SimuSGammaData(inLink,modid,n,p)
%function retSD = SimuSGammaData(inLink,modid,n,p)
%
% Simulate data from a sparse gamma model.
%
% Model: 1/E[Y] = mu + X*beta
%
% Input:
%   inLink: link function of the generalized linear model;
%   modid: index of the model;
%   n: number of observation in each dataset;
%   p: number of variables in X.
%
% Output:
%   retSD: = [Y,X,mu,beta,p]
%
% Dabao Zhang     November 4, 2011
%
% Revised by:
%
% Copyright (c) 2011-12 by Dabao Zhang.

error('This gamma model is NOT implemented!');

retSD.X = zeros(n,p);
retSD.Y = zeros(n,1);
retSD.mu = 0;
retSD.beta = zeros(p,1);



%-------------------------------------------------------------------------
function retSD = SimuSInvGaussianData(inLink,modid,n,p)
%function retSD = SimuSInvGaussianData(inLink,modid,n,p)
%
% Simulate data from a sparse inverse Gaussian model.
%
% Model: 1/(E[Y])^2 = mu + X*beta
%
% Input:
%   inLink: link function of the generalized linear model;
%   modid: index of the model;
%   n: number of observation in each dataset;
%   p: number of variables in X.
%
% Output:
%   retSD: = [Y,X,mu,beta,p]
%
% Dabao Zhang     November 4, 2011
%
% Revised by:
%
% Copyright (c) 2011-12 by Dabao Zhang.

error('This inverse Gaussian model is NOT implemented!');

retSD.X = zeros(n,p);
retSD.Y = zeros(n,1);
retSD.mu = 0;
retSD.beta = zeros(p,1);

          
          
%-------------------------------------------------------------------------
function retSD = SimuSQuasiData(inLink,modid,n,p)
%function retSD = SimuSQuasiData(inLink,modid,n,p)
%
% Simulate data from a sparse quasi model.
%
% Model: logit(E[Y]) = mu + X*beta
%
% Input:
%   inLink: link function of the generalized linear model;
%   modid: index of the model;
%   n: number of observation in each dataset;
%   p: number of variables in X.
%
% Output:
%   retSD: = [Y,X,mu,beta,p]
%
% Dabao Zhang     November 4, 2011
%
% Revised by:
%
% Copyright (c) 2011-12 by Dabao Zhang.

error('This quasi model is NOT implemented!');

retSD.X = zeros(n,p);
retSD.Y = zeros(n,1);
retSD.mu = 0;
retSD.beta = zeros(p,1);
