%---------------------------------------MFILE LOGISTIC REGRESSION- AS BERNOULLI------------------------------------------
clear all
close all
 disp('Logistic Regression: Bliss Beetle Data')
 lw = 3; 
 set(0, 'DefaultAxesFontSize', 16);
 fs = 15;
 msize = 10;

%Insects were exposed to gaseous carbon disulphide for a period of 5 hours.
%Eight experiments were run with different concentrations of carbon disulphide.
%Bliss (1935). The calculation of the dosage-mortality curve. Annals of Applied Biology.
%
%ln(concentration)    #killed     #exposed
beetle=[...
3.893 6 59;...
3.970 13 60;...
4.042 18 62;...
4.108 28 56;...
4.171 52 63;...
4.230 53 59;...
4.285 61 62;...
4.338 60 60];

%%%%%%%%%%%%%%%%%

m=length(beetle(:,1)) ;
x=[]; y=[];
for j=1:m % represent Binomial observations as Bernoulli
x=[x,beetle(j,1)*ones(1,beetle(j,3))];
y=[y,ones(1,beetle(j,2)),zeros(1,beetle(j,3)-beetle(j,2))];
end
beta0 = [0 0];
[betamle,stderr, ps, dev] = logisticmle(y,x);
%-----------------------------
lp = betamle(1) + betamle(2) * x; phat = exp(lp)./(1+exp(lp));
deviance = 2 * sum( y.* log((y+eps)./phat)  + ...
    (1-y) .* log(   (1-y+eps)./(1-phat)) );
dres = sign(y - phat) .* sqrt(2 * y.* log((y+eps)./phat)  + ...
    (1-y) .* log(   (1-y+eps)./(1-phat)) );
%----------------------------
figure(1)
xx=3.5:0.02:4.5;  
linpart=betamle(1) + betamle(2)*xx;
pp=exp(linpart)./(1+exp(linpart));
plot(xx,pp, '-')
hold on
plot( beetle(:,1), beetle(:,2)./beetle(:,3),'o','markersize',msize)

figure(2)
plot(dres,'*')
