close all
clear all

ImmPeroxNeg=[...
19, 25, 30, 34, 37, 46, 47, 51, 56, 57, 61, 66, 67, 74, 78, 86,...
122, 123, 130, 130, 133, 134, 136, 141, 143, 148, 151, 152,...
153, 154, 156, 162, 164, 165, 182, 189];
CensorIPN=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,...
  1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1];

ImmPeroxPos =[...
22, 23, 38, 42, 73, 77, 89, 115, 144];
CensorIPP=[0,0,0,0,0,0,0,0,1];

%number of observed (non-censored)
k1 = sum(1-CensorIPN)   %16
k2 = sum(1-CensorIPP)   %8

% MLE estimators of rate lambda for 2 samples.
hatlam1 = k1/sum(ImmPeroxNeg)  %0.0042
hatlam2 = k2/sum(ImmPeroxPos)  %0.0128


varlam1 = k1/(sum(ImmPeroxNeg))^2  %0.0042
varlam2 = k2/(sum(ImmPeroxPos))^2  %0.0128

z0975 = norminv(0.975);
[hatlam1 - z0975*sqrt(varlam1), hatlam1 + z0975*sqrt(varlam1)] 

 exp([log(hatlam1) - z0975*sqrt(1/k1) log(hatlam1) + z0975*sqrt(1/k1)])


[reclambdahat1 lamci1] = mle(ImmPeroxNeg, ...
   'distribution','exponential','censoring',CensorIPN)  
   % 237.6250
   % 153.6769  415.7289
[reclambdahat2 lamci2] = mle(ImmPeroxPos, ...
   'distribution','exponential','censoring',CensorIPP) 
   % 77.8750
   % 43.1959   180.3793
%(MATLAB parametrization) scale to rate
lambdahat1 = 1/reclambdahat1  %0.0042
lambdahat2 = 1/reclambdahat2  %0.0128

z = (log(hatlam1) - log(hatlam2))/sqrt(1/k1 + 1/k2)  %-2.5763
p = normcdf(z)  %0.0050

z0975 = norminv(0.975);
[hatlam1 - z0975*sqrt(varlam1), hatlam1 + z0975*sqrt(varlam1)] 
  %0.0021    0.0063
  
[hatlam1 - z0975* hatlam1/sqrt(k1), hatlam1 + z0975* hatlam1/sqrt(k1)] 
  %0.0021    0.0063
  
  
 exp([log(lambdahat1) - z0975*sqrt(1/k1) ...
      log(lambdahat1) + z0975*sqrt(1/k1)])  
  %0.0026    0.0069
 1./lamci1 %0.0024 0.0065
   


