function retW = GEBCwMLE(inX,inPatt)
%function retW = GEBCwMLE(inX,inPatt)
% Calculate the marginal MLE of the prior weights
%
% The marginal log-likelihood function is defined as folows:
%   LogL(wn,wp)=sum_{i=1}^{p} log[(1-wp-wn)\phi(y_i)+wp*g+(y_i)+wn*g-(y_i)];
%   retW = (wn,wp).
%
% Using Newton-Raphson to solve for wp and wn.
%       wp+3*wn<=2;
%       3*wp+wn<=2;
%
% Note:
%       inPatt=(-1,1,1) implies that wp, wn \in [0,1];
%       inPatt=(1,1,1) implies that wp=wn \in [0,0.5];
%       inPatt=(1,0,1) implies that wp+wn=1, i.e., wp=wn=0.5;
%       inPatt=(0,1,1) implies that wn=0, wp \in [0,2/3];
%       inPatt=(-1,1,0) implies that wp=0, wn \in [0, 2/3];
%
% Dabao Zhang     July 1, 2006
%
% Revised by:
%
% Copyright (c) 2006 by Zhang & Zhang.


if( (inPatt(1)==-1)&(inPatt(2)==1)&(inPatt(3)==1) )
    %inPatt=(-1,1,1) implies that wp, wn \in [0,1]
    [wp,wn] = GEBCUncMLE(inX);    
elseif( (inPatt(1)==0)&(inPatt(2)==1)&(inPatt(3)==1) )
    %inPatt=(0,1,1) implies that wn=0, wp \in [0,2/3];
    wn = 0;
    wp = GEBCPSigMLE(inX);
elseif( (inPatt(1)==-1)&(inPatt(2)==1)&(inPatt(3)==0) )
    %inPatt=(-1,1,0) implies that wn \in [0,2/3], wp=0;
    wp = 0;
    wn = GEBCPSigMLE(-inX);
elseif( (inPatt(1)==1)&(inPatt(2)==1)&(inPatt(3)==1) )
    %inPatt=(1,1,1) implies that wp=wn \in [0,0.5];
    wp = EBCwMLE(inX)/2;   wn = wp;
elseif( (inPatt(1)==1)&(inPatt(2)==0)&(inPatt(3)==1) )
    %inPatt=(-1,0,1) implies that wp+wn=1, i.e., wp=wn=0.5;
    wp = 0.5;   wn = 0.5;
else
    error('No such estimate exists!');
    wp = 0.5;   wn = 0.5;
end

%disp(['wp = ' num2str(wp)]);
%disp(['wn = ' num2str(wn)]);

retW = [wn,wp];



%--------------------------------------------------------------------------
function [wpRet, wnRet] = GEBCUncMLE(inData)
%function [wpRet, wnRet] = GEBCUncMLE(inData)
% Calculate the marginal MLE of the prior weights
%
% The marginal log-likelihood function is defined as folows:
%   LogL(wp, wn)=sum_{i=1}^{p} log[(1-wp-wn)\phi(y_i)+wp*g+(y_i)+wn*g-(y_i)];
%
% Using the MatLab function fmincon to find the optimal wp and wn,
%       wp+3*wn<=2;
%       3*wp+wn<=2;
%
% Dabao Zhang     July 3, 2006
%
% Revised by:
%
% Copyright (c) 2006 by Zhang & Zhang.

A = [3 1; 1 3; -1 0; 0 -1];
b = [2; 2; 0; 0];
tmpW = [1/3;1/3];
tmpOpts = optimset('LargeScale','off','Display','off');
tmpW = fmincon(@(x) GEBCNMLikeL(x,inData),tmpW,A,b,[],[],[],[],[],tmpOpts);
wnRet = tmpW(1);
wpRet = tmpW(2);



%--------------------------------------------------------------------------
function retNMLL = GEBCNMLikeL(inW,inX)
%function retNMLL = GEBCNMLikeL(inW,inX)
% Calculate negative value of the marginal likelihood function
%
% The marginal log-likelihood function is defined as folows:
%   LogL(wp, wn)=sum_{i=1}^{p} log[(1-wp-wn)\phi(y_i)+wp*g+(y_i)+wn*g-(y_i)];
%
% Dabao Zhang     July 3, 2006
%
% Revised by:
%   Dabao Zhang, Feb. 19, 2008 -- Avoid non-real function values
%
% Copyright (c) 2006 by Zhang & Zhang.

% 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));

if( any(inW<0)|(sum(inW)>1) )
    retNMLL = realmax;
else
    %retNMLL = -sum(log((1-inW(1)-inW(2))*tmpDPhi+inW(1)*tmpGn+inW(2)*tmpGp));
    tmpLikeL = (1-inW(1)-inW(2))*tmpDPhi+inW(1)*tmpGn+inW(2)*tmpGp;
    tmpLikeL(tmpLikeL<=eps) = eps;
    retNMLL = -sum(log(tmpLikeL));
end

if(~isreal(retNMLL))
    retNMLL = realmax;
end



%--------------------------------------------------------------------------
function wpRet = GEBCPSigMLE(inData)
%function wpRet = GEBCPSigMLE(inData)
% Calculate the marginal MLE of the prior weight wp (wn=0)
%
% The marginal log-likelihood function is defined as folows:
%   LogL(wp)=sum_{i=1}^{p} log[(1-wp)\phi(y_i)+wp*g+(y_i)];
%
% Using Bisection-Search to find optimal wp (wn=0).
%       wp<=2/3;
%
% Dabao Zhang     Feb. 14, 2005
%
% Revised by:
%       Dabao Zhang, July 1, 2006 -- Speed up the program.
%
% Copyright (c) 2004-2006 by Zhang & Zhang.


% Redefine variable names and initialize variables for easy programming
p = size(inData,1);
tmpGp = zeros(p,1);
nMaxIter = 5000; % maximum number of iterations.
eLimit = 1e-12; % the criteria for convergence.

% specify initial values
wpHi = 2/3;
wpLo = 0;
wpRet = 1/3;

% compute the G function;
tmpCPhi = normcdf(inData);
tmpDPhi = normpdf(inData);
tmpGp = (tmpCPhi-tmpDPhi.*(sqrt(pi/2)+inData))./(inData.*inData*sqrt(pi/2));

% Iteration to searach roots;
curIter = 0;
diff = 1;

while( (diff>=eLimit)&(curIter<nMaxIter) )
    curIter = curIter+1;

    fWpRet = sum((tmpGp-tmpDPhi)./((1-wpRet)*tmpDPhi+wpRet*tmpGp));
    if( fWpRet > 0 )
        wpLo = wpRet;
    elseif( fWpRet < 0 )
        wpHi = wpRet;
    else    % fWpRet == 0
        wpLo = wpRet;
        wpHi = wpRet;
    end
        
    %check convergence;
    diff = abs(wpHi-wpLo);
    wpRet = (wpLo+wpHi)/2;
end

if( curIter>=nMaxIter )
    disp('Searching for marginal MLE of wp does not converge!');
end



%--------------------------------------------------------------------------
function wRet = EBCwMLE(inData)
%function wRet = EBCwMLE(inData)
% Calculate the marginal MLE of the prior weight.
%
% The marginal log-likelihood function is defined as folows:
%   LogL(w)=sum_{i=1}^{p} log[(1-w)\phi(y_i) + w*g(y_i)];
%
% Using Newton-Raphson to solve for w.
%
% MinZhang      September 15, 2004
%
% Revised by:
%       Dabao Zhang, July 1, 2006 -- Speed up the program.
%
% Copyright (c) 2004-2006 by Zhang & Zhang.

% Redefine variable names and initialize variables for easy programming
p = size(inData,1);
nMaxIter = 250; % maximum number of iterations.
eLimit = 1e-12; % the criteria for convergence.

% Specify initial values
wRet = (sum(inData>0)+sum(inData<0))/p;
wTmp = wRet;

% Compute the G functions;
tmpDPhi = normpdf(inData);
tmpG = (1/sqrt(2*pi)-tmpDPhi)./(inData.*inData);

% Iteration to searach roots;
diff = 1;
curIter = 0;
while( (diff>=eLimit)&(curIter<nMaxIter) )
    curIter = curIter+1;
    tmpC = (1-wRet)*tmpDPhi+wRet*tmpG;
    fW = sum((tmpG-tmpDPhi)./tmpC);
    dfW = sum(((tmpG-tmpDPhi)./tmpC).^2);
    wTmp = wRet+fW/dfW;
    
    if( wTmp<0 )
        wTmp = wRet/2;
    elseif( wTmp>1 )
        wTmp = 1-(1-wRet)/2;
    end

    % Check convergence;
    diff = abs(wTmp-wRet);
    wRet = wTmp;
end

if( curIter>=nMaxIter )
    disp('Searching marginal MLE for w does not converge!');
end
