function retPV = SimuStudySIGPOCRE(inFamily,inLink,modid,nd,n,p,ns,ptype)
%function retRes = SimuStudySIGPOCRE(inFamily,inLink,modid,nd,n,p,ns,ptype)
%
% Simulation study of the combination of information criterion and pooled
% p-values when gPOCRE is applied to generalized linear models.
%
% 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;
%   nd: number of datasets;
%   n: number of observation in each dataset;
%   p: number of variables in X;
%   ns: number of random splits;
%   ptype -- type of penalization:
%            'ebt' ~ EBT w/o z-transformation,
%            'ebtz' ~ EBT w/ z-transformation,
%            'l1' ~ L1,
%            'scad' ~ SCAD,
%            'mcp' ~ MCP.
%
% Output:
%   retPV: p-values for different model selection criteria, p x nd+3 x 13.
%
% Notes on Family & Link
%   Binomial: logit (default), probit, cloglog
%   Poisson: log (default), identity, sqrt
%   Gaussian: identity (default), log, inverse
%   Gamma: inverse (default), identity, log
%   Inverse Gaussian (IGaussian): 1/mu^2 (default), inverse, identity, log
%   Quasi: logit (default), probit, cloglog, identity, inverse, log,
%          1/mu^2, sqrt
%
% Note: quasi family should be able to specify its own power link function
%
% Dabao Zhang    May 16, 2013
%
% Revised by:
%
% Copyright (c) 2013 by Dabao Zhang.

% nohup matlab -nodisplay -r "SimuStudySIGPOCRE('bin','logit',1,100,100,1000,50,'ebtz');exit" >/dev/null 2>error.log &

error(nargchk(3, 8, nargin));
if( nargin<8 ) ptype = 'ebtz'; end
if( nargin<7 ) ns = 50; end
if( nargin<6 ) p = 1000; end
if( nargin<5 ) n = 100; end
if( nargin<4 ) nd = 100; end

% Initialization
inOption.tol = 1e-6;
inOption.maxit = 100;
inOption.bctype = 1;
inOption.maxcmp = 50;
inOption.maxvar = n/2;
inOption.ptype = ptype;
retPV = ones(p,nd+3,13);    % p x (nd+3) x 13

resFName = [upper(inFamily(1:3)) inLink 'M' num2str(modid) ...
            'N' num2str(n) 'P' num2str(p) '.xls'];

nH = round(n)/2;
msc = [2 3 0:0.1:1];
for j = 1:nd
    disp(['Model' num2str(modid) ': dataset=' num2str(j)]);
    
    rand('state',j*100);
    randn('state',j*100);
    
    % Generate training data
    tmpYX = SimuSGLMData(inFamily,inLink,modid,n,p);
    %dataFName = [upper(inFamily) inLink 'M' num2str(modid) 'N' ...
    %             num2str(n) 'P' num2str(p) 'DS' num2str(j) '.txt'];
    %YX = [tmpYX.Y tmpYX.X];  save(dataFName,'YX','-ascii'); clear YX;
    
    pv = ones(p,ns,13);
    for k=1:ns
        disp(['   Random Split: ' num2str(j)]);

        tmpIdx = randperm(n);
        iIdx = tmpIdx(1:nH);
        oIdx = tmpIdx((nH+1):end);

        % Fit gPOCRE
        tmpRes = gPOCREPath(inY(iIdx,:),inX(iIdx,:),inFamily,inLink,inOption);

        % Find optimal model for different model selection criteria
        for m = 1:13
            tmpRes = SelectModel(tmpRes, msc(m));
            xIdx = find(abs(tmpRes.beta)>eps);

            % Calculate p-values
            nx = length(xIdx);
            if( nx>0 )
                [b,d,ss] = glmfit(inX(oIdx,xIdx),inY(oIdx,:),...
                                  inFamily,'link',inLink);
                pv(xIdx,k,m) = min(ss.p(2:end)*nx,ones(nx,1));
            end
        end        
    end

    % Pool p-values
    for m=1:13
        qpv = zeros(p,19);
        for k = 1:19
            gamma = k*0.05;
            gpv(:,k) = min(quantile(pv(:,j,m),gamma,2)/gamma,ones(p,1));
        end
        retPV(:,j,m) = min((1-log(0.05))*min(gpv,[],2),ones(p,1));
    end    
end

for m=1:13
    retPV(:,nd+1,m) = mean(retPV(:,1:nd,m),2);
    retPV(:,nd+2,m) = median(retPV(:,1:nd,m),2);
    retPV(:,nd+3,m) = std(retPV(:,1:nd,m)')';
end

xlswrite(resFName,retPV(:,:,1),'AIC');
xlswrite(resFName,retPV(:,:,2),'AICc');
xlswrite(resFName,retPV(:,:,3),'BIC');
for k=4:13
    xlswrite(resFName,retPV(:,:,k),['EBIC:' num2str(msc(k))]);
end
disp(['    Congratulations! The job is done!');
