function Leif=EfficientELlikelihood(W,delta,theta);

%%%% This code is made by Dr. W Liang, wliang@xmu.edu.cn

%%%% This function is used to calculate the Efficient-EL likehood
%%%% the estimating equation for this likelihood is E(X-theta)=0. 
%%%%
%%%% input: doubly censored data = (W,delta)
%%%%        parameter mean = theta
%%%% output: Efficient-EL likehood = Leif


%%%% sort W %%%%
[W,I]=sort(W);
delta=delta(I);
delta(n)=1;
delta(1)=1;
%%%% Using the doubly censoed data(W, delta) to solve the Turnbull estimator hatF
n=size(W,1);  
theta1=1-unifcdf(W,0,3);
theta1(n)=0;
eps=1e-8;
Wrow=repmat(W',n,1);     %%% [W';W';W']
Wcolumn=repmat(W,1,n);   %%% [W,W,W]
compareW1=(Wrow>Wcolumn);
compareW2=(Wrow<=Wcolumn);
compareW3=compareW1+diag(ones(n,1));
compareW4=compareW2-diag(ones(n,1));
kk2=1;
theta2=theta1;
dist=1;
Del=-1;
while Del<0;
    theta1=theta2;
    kk1=1;
    while dist>eps;
        theta2 = sum(compareW1,2)/n + theta1.*(compareW2*((delta==2)./(theta1+eps)))/n - (1-theta1).*(compareW1*((delta==3)./(1-theta1+eps)))/n;
        dist=max(theta1-theta2);
        theta1=theta2;
        kk1=kk1+1;
        if kk1>1e5;
            fprintf('Maximum number of iterations reached.\n');
            break;
        end
    end    
    %%%% check condition (2) of Theorem2  %%%%
    Delta=1-mean(compareW4.*(delta==2)'./(theta1+eps)'+compareW3.*(delta==3)'./(1-theta1+eps)',2);
    Del=min(Delta);
    theta2=unifrnd(0,1,n,1);
    theta2=sort(theta2,'descend');
    kk2=kk2+1;
    if kk2>1e5;
        fprintf('Maximum number of iterations reached.\n');
        break;
    end
end
%%%%%  the Turnbull estimator hatF
F=1-theta1;
F(n)=1;
%%%%% the self-consistent estimators hatG_L and hatG_R
G_L=F;
G_R=F;
kfirst=find(F>0,1,'first');
for k=kfirst:1:n;
    G_R(k)=1-mean((W>W(k)).*(delta==3)./(F+eps));
end
if kfirst>1;
    G_R(1:kfirst-1)=0;
end
klast=find(F<1,1,'last');
for k=1:1:klast;   
    G_L(k)=mean((W<=W(k)).*(delta==2)./(1-F+eps));
end
if klast<n;
    G_L(klast+1:n)=1;
end
        
%%%%  Solve the equation (7) and get hat_h   %%%%%
hat_theta=(W-theta); 
for k=kfirst:1:n;
    K2(k)=mean((W>=W(k)).*(delta==3)./(F+eps).^2);
end
if kfirst>1;
    K2(1:kfirst-1)=K2(kfirst);
end
for k=1:1:klast;
    K1(k)=mean((W<W(k)).*(delta==2)./(1-F+eps).^2);
end
if klast<n;
    K1(klast+1:n)=K1(klast);
end

for ik=1:1:n;
    for jk=1:1:n;
        if ik<jk;
            K(ik,jk)=(K1(ik)+K2(jk))/(G_R(jk-1)-G_L(jk-1)+eps)*(delta(jk)==1)/n;
        elseif jk==ik;
            if ik==1;
                K(1,1)=(G_R(1)-0)+(K1(1)+K2(1))/(G_R(1)-0+eps)*(delta(1)==1)/n;
            else
                K(ik,ik)=G_R(ik-1)-G_L(ik-1)+(K1(ik)+K2(ik))/(G_R(ik-1)-G_L(ik-1)+eps)*(delta(ik)==1)/n;
            end
        else
            if jk==1;
                K(ik,1)=(K1(1)+K2(ik))/(G_R(1)-0+eps)*(delta(jk)==1)/n;
            else
                K(ik,jk)=(K1(jk)+K2(ik))/(G_R(jk-1)-G_L(jk-1)+eps)*(delta(jk)==1)/n;
            end
        end
    end
end
tempK=zeros(n,1)+eps;
invK=K+diag(tempK);
hat_h=inv(invK)*hat_theta; 

%%%%  get hat_psi   %%%%%
hat_psi=hat_h;
GZtemp=[G_R(1);G_R(1:n-1)];
GYtemp=[0;G_L(1:n-1)];
for k=kfirst:1:n;
    hat_psi(k)=(delta(k)==1)*hat_h(k)+(delta(k)==2)*(mean((W>W(k)).*(delta==1).*hat_h./(GZtemp-GYtemp+eps)))/(1-F(k)+eps)+(delta(k)==3)*(mean((W<=W(k)).*(delta==1).*hat_h./(GZtemp-GYtemp+eps)))/(F(k)+eps);
end
if kfirst>1;
    hat_psi(1:kfirst-1)=hat_psi(kfirst);
end
if klast<n;
    hat_psi(klast+1:n)=hat_psi(klast);
end

%%%% Using standard EL function to get the efficient-EL likelihood %%%%%
[logelr, lambda, grad, hess, wts, nits]=elm(hat_psi);
Leif=-2*logelr;

 