function retSD = SimuSGLMData(inFamily,inLink,modid,n,p,inTheta)
%function retSD = SimuSGLMData(inFamily,inLink,modid,n,p,inTheta)
%
% 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;
%   inTheta: scale of coefficients (=1 by default).
%
% Output:
%   retSD = [Y,X,mu,beta]
%
% Dabao Zhang     November 7, 2011
%
% Revised by:
%       Dabao Zhang, May 21, 2013 -- Family & Link.
%       Dabao Zhang, Feb. 12, 2016 --- Use inTheta for BC models.
%
% Copyright (c) 2011-12 by Dabao Zhang.

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

if( nargin<6 ) inTheta = 1; end

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

switch inFamily
    case 'binomial'
        retSD = SimuSBinomialData(inLink,modid,n,p);
    case 'multibc'
        retSD=SimuSMultiBCData(inLink,modid,n,p,inTheta);
    case 'multipo'
        retSD=SimuSMultiPOData(inLink,modid,n,p);
    case 'multicr'
        retSD=SimuSMultiCRData(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);
    case 'cox'
        retSD = SimuSCoxData(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;
        retSD.beta(1:2) = 1;  retSD.beta(103:104) = -1; retSD.beta(205:206) = 1;
        if( modid==1 ) rho = 0.3; else rho = 0.6; 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*retSD.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*retSD.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 = retSD.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*retSD.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);


%-------------------------------------------------------------------------
function retSD = SimuSMultiBCData(inLink,modid,n,p,inTheta)
%function retSD = SimuSMultiBCData(inLink,modid,n,p,inTheta)
%
% Simulate data from a sparse multinomial baseline categorical model.
%
% Model: multilogit(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;
%   inTheta: scale of the coefficients (=1 by default).
%
% Output:
%   retSD: = [Y,X,mu,beta,p,q]
%
% Min Ren     March 17, 2015
%
% Revised by:
%       Dabao Zhang, Feb. 12, 2016 --- Take inTheta.
%
% Copyright (c) 2011-12 by Dabao Zhang.

if( nargin<5 ) inTheta = 1; end

% q is the number of the categories minus 1
q = 2;

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

switch modid
    case {1,2}
        delta = 1;

        %retSD.beta(1:10,:) = 2;  retSD.beta(101:110,:) = 1;
        retSD.beta(1:5,:)=2*delta;
        retSD.beta(6:10,:)=-2*delta;
         %retSD.beta(1:2,:)=100;
         %retSD.beta(3:4,:)=-50;
        retSD.beta(101:105,1)=2*delta;
        retSD.beta(201:205,2)=2*delta;
        if( modid==1 ) rho = 0.2; 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);

    case 5 % with only common components
        retSD.beta(1:10,1) = inTheta;
        retSD.beta(101:110,1) = inTheta;
        retSD.beta(1:10,2) = .5*(1-sqrt(3))*inTheta;
        retSD.beta(101:110,2) = .5*(1+sqrt(3))*inTheta;

        rho = 0.5;
        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 6 % with only common components
        retSD.beta(1:10,1) = sqrt(3)*inTheta;
        retSD.beta(101:110,1) = inTheta;
        retSD.beta(1:10,2) = sqrt(3)*inTheta;
        retSD.beta(101:110,2) = -inTheta;

        rho = 0.5;
        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 7 % Need to compare the power of detecting common vs. peculiar components
        retSD.beta(1:10,1) = inTheta;
        retSD.beta(1:10,2) = inTheta;
        retSD.beta(101:110,1) = inTheta;
        retSD.beta(201:210,2) = inTheta;

        rho = 0.5;
        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;
        
    otherwise
        error('This multinomial model is NOT implemented!');
end

prob = arrayfun(@(i)(ginv(eta(i,:),inLink)),[1:n]','UniformOutput',false);
prob = cell2mat(prob);

tmpSum=sum(prob,2);
P=horzcat(1-tmpSum,prob);
P(P<0)=0;
P = bsxfun(@rdivide,P,sum(P,2));

retSD.Y = mnrnd(ones(n,1),P);


%-------------------------------------------------------------------------
function retSD = SimuSMultiPOData(inLink,modid,n,p)
%function retSD = SimuSMultiPOData(inLink,modid,n,p)
%
% Simulate data from a sparse multinomial proptional odds model.
%
% Model: pomultilogit(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,q]
%
% Min Ren     March 17, 2015
%
% Revised by:
%
% Copyright (c) 2011-12 by Dabao Zhang.

% q is the number of the categories minus 1
q = 2;

retSD.X = zeros(n,p);
retSD.Y = zeros(n,q);
%retSD.mu = cumsum(normrnd(0,1,[1 q]));% generate the intercept
retSD.mu = [-log(8),log(9)];% for q=2
retSD.beta = zeros(p,1);

switch modid
    case {1,2}
        %retSD.beta(1:10) = 1;  retSD.beta(101:110) = -1;
        retSD.beta(1:5)=-1.5;retSD.beta(6:10)=-.5;
        if( modid==1 ) rho = 0.2; 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 = ones(n,1)*retSD.mu+retSD.X*retSD.beta*ones(1,q);
        
    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 = ones(q,1)*retSD.mu+retSD.X*retSD.beta*ones(1,q);
        
    otherwise
        error('This multinomial model is NOT implemented!');
end

prob = arrayfun(@(i)(ginv(eta(i,:),inLink)),[1:n]','UniformOutput',false);
prob = cell2mat(prob);

tmpSum=sum(prob,2);
P=horzcat(1-tmpSum,prob);
P(P<0)=0;
P = bsxfun(@rdivide,P,sum(P,2));

retSD.Y = mnrnd(ones(n,1),P);


%-------------------------------------------------------------------------
function retSD = SimuSMultiCRData(inLink,modid,n,p)
%function retSD = SimuSMultiCRData(inLink,modid,n,p)
%
% Simulate data from a sparse multinomial continuation ratio model.
%
% Model: pomultilogit(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,q]
%
% Min Ren     March 19, 2015
%
% Revised by:
%
% Copyright (c) 2011-12 by Dabao Zhang.

% q is the number of the categories minus 1
q = 2;

retSD.X = zeros(n,p);
retSD.Y = zeros(n,q);
%retSD.mu = cumsum(normrnd(0,1,[1 q]));% generate the intercept
retSD.mu = [-log(8),log(9)];% for q=2
retSD.beta = zeros(p,1);

switch modid
    case {1,2}
        %retSD.beta(1:10) = 1;  retSD.beta(101:110) = -1;
        retSD.beta(1:5)=-2;retSD.beta(6:10)=1;
        if( modid==1 ) rho = 0.2; 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 = ones(n,1)*retSD.mu+retSD.X*retSD.beta*ones(1,q);
        
    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 = ones(q,1)*retSD.mu+retSD.X*retSD.beta*ones(1,q);
        
    otherwise
        error('This multimial model is NOT implemented!');
end

prob = arrayfun(@(i)(ginv(eta(i,:),inLink)),[1:n]','UniformOutput',false);
prob = cell2mat(prob);

tmpSum=sum(prob,2);
P=horzcat(1-tmpSum,prob);
P(P<0)=0;
P = bsxfun(@rdivide,P,sum(P,2));

retSD.Y = mnrnd(ones(n,1),P);


%-------------------------------------------------------------------------
function retSD = SimuSCoxData(modid,n,p)
%function retSD = SimuSCoxData(modid,n,p)
%
% Simulate data from a sparse proportional hazards model.
%
% Model: log(h(t)/h0(t)) = X*beta
%
% Input:
%   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,q]
%
% Dabao Zhang     April 4, 2015
%
% Revised by:
%
% Copyright (c) 2011-12 by Dabao Zhang.

retSD.X = zeros(n,p);
retSD.Y = zeros(n,1);
retSD.mu = 0; % always zero in Cox Proportional Hazards Model
retSD.beta = zeros(p,1);

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

switch modid
    case {1,2} % High/Mild Correlation
        retSD.beta(1:110) = [ones(10,1); zeros(90,1); 0.5*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
        eta = retSD.mu+retSD.X*retSD.beta;
        
    case 3 % Clusering
        retSD.beta(1:30) = .75;

        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.X*retSD.beta-5;
        
    case 4 % Measurement Error Model
        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 = 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 model is NOT supported!');
end

expEta = exp(eta);
retSD.Y = wblrnd(expEta,ones(n,1)*3);
tmpT = quantile(retSD.Y(:,1),0.95);
retSD.Y(:,2) = (retSD.Y(:,1)>tmpT);
retSD.Y(retSD.Y(:,1)>tmpT,1) = tmpT;
