function [f,g] = LL_poiss(Y, Xa,Xnb, EstimOpt,b)

NB = EstimOpt.NB;
Censored = EstimOpt.Censored;
NVarA = EstimOpt.NVarA;
NVarNB = EstimOpt.NVarNB;
RealMin = EstimOpt.RealMin; 

if NB == 0
    if Censored == 0
        f = (Xa*b).*Y - gammaln(Y+1) - exp(Xa*b);
    else
        f = zeros(length(Y),1);
        f(Y ~= Censored) = (Xa(Y ~= Censored,:)*b).*Y(Y ~= Censored) - gammaln(Y(Y ~= Censored)+1) - exp(Xa(Y ~= Censored,:)*b);
        cXa = Xa(Y == Censored,:);
        lam = exp(cXa*b);
        vec = 0:(Censored - 1);
        vec = vec(ones(length(lam),1),:);
        fitc = cXa*b;
        fitc = fitc(:,ones(Censored,1));
     
        %ftmp = 1 - exp(-lam).*sum(exp(vec.*fitc - gammaln(vec+1)),2);
        ftmp = 1 - sum(exp(vec.*fitc - lam - gammaln(vec+1)),2);
       % find(isnan(ftmp))
     %   any(ftmp == 0)
        %ftmp = 1 - sum(exp(vec.*fitc - gammaln(vec+1)-lam(:,ones(Censored,1))),2);
        if RealMin == 0
            f(Y == Censored) = log(ftmp);
        else
            ftmp = max(ftmp,realmin);
            f(Y == Censored) = log(ftmp);
        end
    end
elseif NB == 1
    theta = exp(-Xnb*b(NVarA+1:NVarA+NVarNB));
    fit = Xa*b(1:NVarA);
    if Censored == 0
        f = gammaln(theta+Y)-gammaln(theta)-gammaln(Y+1)-theta.*(Xnb*b(NVarA+1:NVarA+NVarNB))+Y.*fit - (Y+theta).*log(theta+exp(fit));
    else
        f = zeros(length(Y),1);
        f(Y ~= Censored) = gammaln(theta(Y ~= Censored)+Y(Y ~= Censored))-gammaln(theta(Y ~= Censored))-gammaln(Y(Y ~= Censored)+1)...
            -theta(Y ~= Censored).*(Xnb(Y ~= Censored,:)*b(NVarA+1:NVarA+NVarNB))+Y(Y ~= Censored).*fit(Y ~= Censored) - (Y(Y ~= Censored)+theta(Y ~= Censored)).*log(theta(Y ~= Censored)+exp(fit(Y ~= Censored)));

        cfit = fit(Y == Censored);
        ctheta = theta(Y == Censored);
        u = ctheta./(ctheta+exp(cfit));
        ut = (u.^ctheta)./gamma(ctheta);
        
        vec = 0:(Censored - 1);
        vec = vec(ones(length(cfit),1),:);
        cfit = cfit(:,ones(Censored,1));
        ctheta = ctheta(:,ones(Censored,1));
        ftmp = sum(exp(gammaln(ctheta+vec)-gammaln(vec+1)+vec.*cfit - vec.*log(ctheta+exp(cfit))),2);
        
        ftmp = 1 - ftmp.*ut;
        if RealMin == 0
            f(Y == Censored) = log(ftmp);
        else
            ftmp = max(ftmp,realmin);
            f(Y == Censored) = log(ftmp);
        end
%         if theta + Censored -1< 171
%             u = theta./(theta+lam);
%             u2 = u(:, ones(Censored,1));
%             f(Y == Censored) = log(1 - ((u.^theta)./gamma(theta)).*sum((gamma(theta+vec)./gamma(vec+1)).*(1-u2).^vec,2));
%         else % if Theta + censoring is too large the gamma function became INF and cause problems. We use the fact that for large theta NB2 becomes approximetly Poisson 
%             lam2 = lam(:, ones(Censored,1));
%             f(Y == Censored) = log(1 - exp(-lam).*sum((lam2.^vec)./gamma(vec+1),2));
%         end
    end
elseif NB == 2
    theta = exp(-Xnb*b(NVarA+1:NVarA+NVarNB));
    Q = 2- b(NVarA+NVarNB+1);
    fit = Xa*b(1:NVarA);
    theta = theta.*exp(Q*fit);
    fx = gammaln(theta+Y)-gammaln(theta)-gammaln(Y+1)-theta.*(Xnb*b(NVarA+1:NVarA+NVarNB))+theta*Q.*fit+Y.*fit - (Y+theta).*log(theta+exp(fit));
    if Censored == 0
        f = fx;
    else
        f = zeros(length(Y),1);
        f(Y ~= Censored) = fx(Y ~= Censored);
        
        cfit = fit(Y == Censored);
        ctheta = theta(Y == Censored);
        cXnb = Xnb(Y == Censored,:);

        vec = 0:(Censored - 1);
        vec = vec(ones(length(cfit),1),:);
        cfit = cfit(:,ones(Censored,1));
        ctheta = ctheta(:,ones(Censored,1));
        
        ftmp = sum(exp(gammaln(ctheta+vec)-gammaln(ctheta)-gammaln(vec+1)-ctheta.*(cXnb*b(NVarA+1:NVarA+NVarNB))+ctheta*Q.*cfit+vec.*cfit - (vec+ctheta).*log(ctheta+exp(cfit))),2);
        ftmp = 1 - ftmp;
        if RealMin == 0
            f(Y == Censored) = log(ftmp);
        else
            ftmp = max(ftmp,realmin);
            f(Y == Censored) = log(ftmp);
        end
        
    end
end

f = -f;

if nargout == 2
    Lam = exp(Xa*b(1:NVarA));
    if NB == 0
        if Censored == 0
            g = Xa.*Y(:, ones(1, NVarA)) - Lam(:, ones(1, NVarA)).*Xa;
            g = -g;
        else
            g = zeros(length(Y),EstimOpt.NVarA);
            g(Y ~= Censored,:) = Xa(Y ~= Censored,:).*Y(Y ~= Censored, ones(1, NVarA)) - Lam(Y ~= Censored, ones(1, NVarA)).*Xa(Y ~= Censored,:);
            
            %gtmp = sum(exp(vec.*fitc-gammaln(vec+1)).*(lam(:,ones(Censored,1))-vec),2);
            %gtmp = exp(-lam).*gtmp./ftmp;
            gtmp = sum(exp(vec.*fitc -lam-gammaln(vec+1)).*(lam(:,ones(Censored,1))-vec),2);
            gtmp = gtmp./ftmp;
            
            %gtmp = sign(gtmp).*exp(log(abs(gtmp))-lam - log(ftmp));
            g(Y == Censored,:) = gtmp(:, ones(1, NVarA)).*cXa;
            g = -g;
        end
    else
        if Censored == 0
            g = Xa.*Y(:, ones(1, NVarA)) - (Y(:, ones(1, NVarA)) + theta(:, ones(1, NVarA))).*Lam(:, ones(1, NVarA)).*Xa./(theta(:, ones(1, NVarA))+Lam(:, ones(1, NVarA))) ;
            txnb = theta(:, ones(1, NVarNB)).*Xnb;
            fit2 = Xnb*b(NVarA+1:NVarA+NVarNB);
            gnb  = -psi(theta+Y)+psi(theta)-1+fit2+log(theta+Lam)+ (theta+Y)./(theta+Lam);
            gnb  = txnb.*gnb(:, ones(1, NVarNB));
        else
            txnb = theta(:, ones(1, NVarNB)).*Xnb;
            fit2 = Xnb*b(NVarA+1:NVarA+NVarNB);
            % without censoring
            g = zeros(length(Y),EstimOpt.NVarA);
            g(Y ~= Censored,:) = Xa(Y ~= Censored,:).*Y(Y ~= Censored, ones(1, NVarA)) - (Y(Y ~= Censored, ones(1, NVarA))...
                + theta(Y ~= Censored, ones(1, NVarA))).*Lam(Y ~= Censored, ones(1, NVarA)).*Xa(Y ~= Censored,:)./(theta(Y ~= Censored, ones(1, NVarA))+Lam(Y ~= Censored, ones(1, NVarA))) ;
            
            gnb = zeros(length(Y),EstimOpt.NVarNB);
            gnb_tmp  = -psi(theta(Y ~= Censored)+Y(Y ~= Censored))+psi(theta(Y ~= Censored))-1+fit2(Y ~= Censored)...
                +log(theta(Y ~= Censored)+Lam(Y ~= Censored))+ (theta(Y ~= Censored)+Y(Y ~= Censored))./(theta(Y ~= Censored)+Lam(Y ~= Censored));
            gnb(Y ~= Censored,:)  = txnb(Y ~= Censored,:).*gnb_tmp(:, ones(1, NVarNB));
            
            % With censoring
            gtmp = sum((exp(cfit)-vec).*exp(gammaln(ctheta+vec)-gammaln(vec+1)+vec.*cfit - vec.*log(ctheta+exp(cfit))),2);
            gtmp = gtmp.*ut.*u./ftmp;
            g(Y == Censored,:) = gtmp(:, ones(1, EstimOpt.NVarA)).*Xa(Y == Censored,:);
            
            
%            gtmp = (-psi(ctheta(:,1)) + log(ctheta(:,1)) + 1 - log(ctheta(:,1)+exp(cfit(:,1)))-u).*sum(exp(gammaln(ctheta+vec)-gammaln(vec+1)+vec.*cfit - vec.*log(ctheta+exp(cfit))),2);
%             gtmp2 = psi(ctheta+vec).*exp(gammaln(ctheta+vec) - gammaln(vec+1)+vec.*cfit - vec.*log(ctheta+exp(cfit)))- exp(gammaln(ctheta+vec)-gammaln(vec+1)+vec.*cfit - (vec+1).*log(ctheta+exp(cfit))).*vec;
%             gtmp = gtmp - sum(gtmp2,2);
            gtmp = exp(gammaln(ctheta+vec) - gammaln(vec+1) +vec.*cfit - vec.*log(ctheta+exp(cfit)));
            gtmp = sum(gtmp.*(psi(ctheta+vec) - psi(ctheta) - vec./(ctheta + exp(cfit)) + log(ctheta)- log(ctheta+exp(cfit)) + 1 - exp(log(ctheta) - log(ctheta+exp(cfit)))),2);
            gtmp = gtmp.*ut.*ctheta(:,1)./ftmp;
            gnb(Y == Censored,:) = gtmp(:, ones(1, EstimOpt.NVarNB)).*Xnb(Y == Censored,:);
        end
        g = -[g, gnb];
    end
   
end
