function retRes = RBootsGPOCRE(inY,inX,inRes,inFamily,inLink,nBoots,inOption)
%function retRes = RBootsGPOCRE(inY,inX,inRes,inFamily,inLink,nBoots,inOption)
%
% Residual bootstrap method is combined with generalized Penalized
% Orthogonal-Component Regression (gPOCRE) to evaluate the variable
% significance. Indeed, used is the the one-step bootstrap method with
% Pearson residuals, see
%
% Moulton LH & Zeger SL (1991) Bootstrapping generalized linear models.
% Computational Statistics & Data Analysis 11: 53-63.
% 
% The tuning parameter may be selected based on (a) information criteria,
% such as (E)BIC, AIC, or AICC; (b) cross-validation; (c) a prespecified
% value. By default, EBIC is used to choose the tuning parameter (note that
% EBIC's own tuning parameter is set as 1-.5*log(n)/log(p-1)).
%
% Input:
%   inY: Y response values of original data;
%   inX: X;
%   inRes: result from POCRE, which includes inRes.mu & inRes.beta;
%   inFamily: the distribution family, such as 'normal', 'binomial',
%             'poisson', 'gamma', 'inverse gaussian', 'quasi' (not
%             supported now), 'multibc', 'multipo','multicr';
%   inLink: the link function g(.), such as 'identity' (default for
%           'normal'), 'log' (default for 'poisson'), 'logit' (default for
%           'binomial'), 'probit','comploglog','reciprocal','loglog'
%           (default for 'gamma'),p (numeric, default for 'inverse
%           gaussian'), user-defined link function (not supported now);
%   nBoots: number of bootstrap data sets;
%   inOption:
%       tol -- tolerance,
%       maxit -- maximum number of iterations,
%       maxcmp: maximum number of components;
%       maxvar -- maximum number of selected variables;
%       ptype -- type of penalization,
%                'ebt' ~ EBT w/o z-transformation,
%                'ebtz' ~ EBT w/ z-transformation,
%                'l1' ~ L1,
%                'scad' ~ SCAD,
%                'mcp' ~ MCP;
%       delta -- pace for choosing tuning parameters.
%
% Output:
%   retRes: mnzbeta,sdnzbeta,nnzbeta,pvalue.
%
% Dabao Zhang    July 27, 2016
%
% Revised by:
%
% Copyright (c) 2016 by Dabao Zhang.

error(nargchk(4, 7, nargin));

if( nargin<7 ) inOption = []; end
if( ~isfield(inOption,'tol') ) inOption.tol = 10e-6; end
if( ~isfield(inOption,'maxit') ) inOption.maxit = 100; end
if( ~isfield(inOption,'maxcmp') ) inOption.maxcmp = 50; end
if( ~isfield(inOption,'maxvar') ) inOption.maxvar = size(inX,1)/2; end
if( ~isfield(inOption,'ptype') ) inOption.ptype = 'ebtz'; end

if( (nargin<6)||isempty(nBoots) ) nBoots = 100; end
if( nargin<5 ) inLink = []; 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 = 'cox';
        otherwise
            error('  RBootsGPOCRE: The family is not supported...');
    end
end

[n,p] = size(inX);    % X with n objects m variables;
k = size(inY,2);    % Y with n objects k variables;
if( n~=size(inY,1) )
    error('  RBootsPOCRE: inX and inY should have the same number of rows!');
end

% Initial values:
retRes.mnzbeta = zeros(p,k);
retRes.sdnzbeta = zeros(p,k);
retRes.nnzbeta = zeros(p,k);

% bootstrap estimates of beta are saved into files:
fname = [num2str(fix(clock))];
fname = [fname(~isspace(fname)) 'beta'];

% Calculate the pearson residuals
eta = bsxfun(@plus,inX*inRes.beta,inRes.mu);
tmpHV = sqrt(CalcIRLSResponseVar(eta,inFamily,inLink));
perr = bsxfun(@times,inY-eta,tmpHV);
tmpEta = bsxfun(@rdivide,eta,tmpHV);
tmpX = bsxfun(@rdivide,inX,tmpHV);

% Bootstrap
for j=1:nBoots
    disp(['   Bootstrap Data Set: ' num2str(j)]);

    tmpY = tmpEta + perr(randsample(n,n,1),:);
    tmpRes = POCREPath(tmpY,tmpX,inOption);
    tmpRes = SelectModel(tmpRes);
    
    for m=1:k
        dlmwrite([fname num2str(m) '.txt'],tmpRes.beta(:,m)','-append');
    end
    retRes.mnzbeta = retRes.mnzbeta+tmpRes.beta/nBoots;
    retRes.sdnzbeta = retRes.sdnzbeta+tmpRes.beta.^2/nBoots;
    retRes.nnzbeta = retRes.nnzbeta+(abs(tmpRes.beta)>eps);
end
retRes.mnzbeta = retRes.mnzbeta.*(nBoots./max(retRes.nnzbeta,1));
retRes.sdnzbeta = retRes.sdnzbeta.*(nBoots./max(retRes.nnzbeta,1));
retRes.sdnzbeta = retRes.sdnzbeta - retRes.mnzbeta.^2;
retRes.sdnzbeta(retRes.sdnzbeta<0) = 0;

% Assume est_beta ~ delta_0 * (1-w) + w* N(mnzbeta,sdnzbeta) to calculate
%   p-value = (1-w)+w*Phi(-|mnzbeta|/sdnzbeta).   
retRes.pvalue = retRes.nnzbeta/nBoots;
retRes.pvalue = 1-retRes.pvalue+retRes.pvalue.*normcdf(0,...
                abs(retRes.mnzbeta),retRes.sdnzbeta);
retRes.pvalue(retRes.sdnzbeta<eps) = 1-1/nBoots;
