function retEst = GEBCPostMeds(inX,inPPn,inPPp)
%function retEst = GEBCPostMeds(inX,inPPn,inPPp)
% Compute the posterior medians for the empirical Bayes shrinkage estimate.
% The Newton-Raphson method is used
%
% inData -- has already been standardized.
%
% Min Zhang     September 12, 2004
%
% Revised by:
%       Dabao Zhang, July. 1, 2006 -- Speed the codes;
%       Min Ren, Mar. 5, 2016 -- Modify for  compiling mex file.
%
% Copyright (c) 2004-2006 by Zhang & Zhang.

% Define variables used in the program;
p = size(inX,1);
nMaxIter = 500; % maximum number of iterations.
eLimit = 1.0e-12; % the criteria for convergence.

retEst = zeros(p,1);
tmpEst = 0;

% compute the G functions;
tmpX2 = inX.*inX;
tmpCPhi = normcdf(inX,0,1);
tmpDPhi = normpdf(inX,0,1);
tmpGp = (tmpCPhi-tmpDPhi.*(sqrt(pi/2)+inX))./(tmpX2*sqrt(pi/2));
tmpGn = (1-tmpCPhi-tmpDPhi.*(sqrt(pi/2)-inX))./(tmpX2*sqrt(pi/2));

% Check whether the postMedian is postive, negative or zero first,
%   then calculate the posterior median correpondingly.

for j=1:p
    bOver = 1;
    diff = 1;
    retEst(j) = inX(j);
    
    yL = 0;
    yU = 0;
    if( inPPp(j)>0.5 )
        % CASE I: check the posibility of being positive values;
        yU = max(0,inX(j));
    elseif( inPPn(j)>0.5 )
        % CASE II: check the possibility of beging negative values;
        yL = min(0,inX(j));
    else
        % CASE III: if none of the above is true, the posterior median is set to be zero;
        retEst(j) = 0;
        diff = 0;
    end
    
    if( ((inX(j)>=7)&(inPPp(j)>0.5))|((inX(j)<=-7)&(inPPn(j)>0.5)) )
        retEst(j) = inX(j)-2/inX(j);
        diff = 0; 
    end
    
    curIter = 0;
    while( (diff>=eLimit)&(curIter<nMaxIter) )
        curIter = curIter+1;
        
        if( inPPp(j)>0.5 )
            % CASE I: check the posibility of being positive values;
            tmpEst = 0.5*(yL+yU);
            bOver = (normcdf(inX(j)-tmpEst)-inX(j)*normpdf(inX(j)-tmpEst)+...
                        +(tmpEst*inX(j)-1)*(1-normcdf(tmpEst))...
                        *exp((tmpEst-inX(j)/2)*inX(j)))/(tmpX2(j)*sqrt(pi/2))...
                        -tmpGp(j)/(2*inPPp(j));
            
            diff = abs(retEst(j)-tmpEst);
            retEst(j) = tmpEst;

            yL = tmpEst + (yL-tmpEst)*(bOver<0);
            yU = tmpEst + (yU-tmpEst)*(bOver>0);
        elseif( inPPn(j)>0.5 )
            % CASE II: check the possibility of beging negative values;
            tmpEst = 0.5*(yL+yU);
            %tmpp = normcdf(tmpEst-inX(j));
            %bOver = 2*tmpp/(tmpX2(j)*sqrt(2*pi))+2*tmpp/(inX(j)*sqrt(2*pi))...
            %            +2*(tmpEst*inX(j)-1)*exp(tmpEst*tmpEst/2)*(1-tmpp)...
            %            *normcdf(tmpEst)/tmpX2(j)-tmpGn(j)/(2*inPPn(j));
            bOver = (normcdf(tmpEst-inX(j))+inX(j)*normpdf(tmpEst-inX(j))...
                        +(tmpEst*inX(j)-1)*normcdf(tmpEst)...
                        *exp((tmpEst-inX(j)/2)*inX(j)))/(tmpX2(j)*sqrt(pi/2))...
                        -tmpGn(j)/(2*inPPn(j));
    
            diff = abs(retEst(j)-tmpEst);
            retEst(j) = tmpEst;

            yL = tmpEst + (yL-tmpEst)*(bOver>0);
            yU = tmpEst + (yU-tmpEst)*(bOver<0);
        else    % Should not come to here!
            disp('Something is wrong here!');
        end
    end
    
    if( curIter>=nMaxIter )
        disp('Searching for Posterior Medians does not converge!');
    end
end
