function retPV = SIGPOCRE(inY,inX,inFamily,inLink,nSplits,inOption)
%function retPV = SIGPOCRE(inY,inX,inFamily,inLink,nSplits,inOption)
%
% Penalized Orthogonal-Component Regression for generalized linear
% models (gPOCRE) with significance inference.  The tuning parameter is
% selected based on a prespecified information criterion, such as (E)BIC,
% AIC, or AICC. Please refer gPOCRE to
%
% Lin, Y., Zhang, M., and Zhang, D. (2011) Penalized orthogonal-component
% regression for high-dimensional generalized linear models. Purdue
% Statistics Technical Report ???.
%
% EBIC=-2*log(L)+p*log(n)+2*p*gamma*log(p_T) can be referred to
%
% Chen J and Chen Z (2010) Extended BIC for small-n-large-p sparse GLM.
% Statistica Sinica. doi: 10.5705/ss.2010.216.
%
% The p-value calculation can be referred to
%
% Meinshausen N, Meier L and Buhlmann P (2009) p-values for
% high-dimensional regression. JASA 104: 1671-1681.
%
% Input:
%   inY: Y response values of original data;
%   inX: X;
%   inFamily: the distribution family;
%   inLink: the link function g(.);
%   nSplits: number of random splits;
%   inOption:
%       tol -- tolerance,
%       maxit -- maximum number of iterations,
%       bctype -- type of bias correction: 0 ~ no bias correction,
%                 1 ~ approximate, type 2 -- optimal;
%       maxcmp: maximum number of components;
%       maxvar -- maximum number of selected variables;
%       W: weight matrix;
%       msc -- model selection criterion,
%              0 ~ BIC,
%              (0,1] ~ coefficient used in EBIC (Chen & Chen, 2010),
%              2 ~ AIC (default),
%              3 ~ AICc;
%       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 value of the tuning parameter
%
% Output:
%   retPV: p-values.
%
% Notes on suggested Family & Link combinations,
%   'normal': identity (default), log, -1
%   'binomial': logit (default), probit, comploglog
%   'poisson': log (default), identity, 0.5
%   'gamma': -1 (default), loglog, identity, log
%   'inverse gaussian': -2 (default), -1, identity, log
%   'quasi': logit (default), probit, comploglog, identity, log, loglog,
%            -1,-2,0.5
%
% Note: quasi family should be able to specify its own power link function
%
% Note: gPOCRE stops the construction of new components when number of
% non-zero-coefficient variables exceeds maxvar (=n/2 by default).
%
% Dabao Zhang    May 15, 2013
%
% Revised by:
%
% Copyright (c) 2013 by Dabao Zhang.

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

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

if( (nargin<5)||isempty(nSplits) ) nSplits = 50; end

if( (nargin<4)||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('  SIGPOCRE: 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(rowY,1) )
    error('  SIGPOCRE: inX and inY should have the same number of rows!');
end

% Only implemented for single response !
if( k>1 )
    error('  gPOCRE4pv: Only implemented for single response!');
end

rand('state',999);
nH = round(n)/2;
pv = ones(p,nSplits);
for j=1:nSplits
    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
    tmpRes = SelectModel(tmpRes, inOption.msc);
    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,j) = min(ss.p(2:end)*nx,ones(nx,1));
    end
end

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