clc; clear all;
antal=importdata('celle_data.txt');
ncoli=antal(:); % antal bakterier i hver firkant som en liste
edges=[4:32]-0.5; % størrelsen af intervallet. 4 er det mindste tal i antal og 32 det største
n=224; % Antal observationer/kasser i alt

figure
    h=histogram(ncoli,edges)
    N=h.Values % hvert tal svarer til højden i histogrammet ved den indgang i søjlen
    % altså hvor mange kasser med et bestemt antal bakterier
    title('Histogram');
    xlabel('Antal colibakterier i en kasse');
    ylabel('Antal kasser');
    hold on
    x=[0:35];
    y=224*poisspdf(x,mean(ncoli));
    plot(x,y)
    hold off
 close all

 x=sum(ncoli); % Antal bakterier i alt - dvs 3324
 t=224; % værdien af t er nok forkert. hjææælp <-----------------------
 lambda_hat=x/t; % Skøn af rateparameter = 14.8393


 p_js = [sum((lambda_hat.^(1:4))./(1:4).*exp(-1*lambda_hat)), ((lambda_hat.^(5:30))./(5:30).*exp(-1*lambda_hat))]
 p_js = [p_js, 1-sum(p_js)]

 ej=224.*p_js; % forventet antal

 % Da jeg ikke kan finde ud af at finde skønnet for rateparameteren ud
 % kan jeg ikke slå kasser sammen, så p.t er der 28 kasser

 %Gx = 2*sum(N.*log(N./ej))
 % den 25 & 26 indgang i N er lig 0, så vi får nogle NaN'er
 % det bliver nok fikset når jeg finder ud af at finde skønnet af
 % rateparameteren og dermed kan ligge kasser sammen

 figure
    title('Histogram');
    xlabel('j');
    ylabel('r_j');
    hold on
    j=[1:length(ej)];
    y=[]
    for i=j
        if i~=25||i~=26
            %y=[y ((N(i)-ej(i))/sqrt(ej(i)))];
        end
    end
    %plot(j,y)
    refline(0,0)
    hold off
close all