% luria0.m
%
% Demonstrates how to do exact Bayesian inference for the Luria-Delbruck problem
%              (c) David J.C. MacKay  Sun 12/1/03
%              This software is distributed under the GNU general public license

clear
more off

data = [ 1,0,3,0,0,5,0,5,0,6,107,0,0,0,1,0,0,64,0,35 ];
% total number of cells
N = 2.4e8 ;

% find largest value of n needed
nmax=max(data); rmax=nmax;

alln = (1:nmax) ;
allzn = (0:nmax) ;

% define assumed distribution of n given one mutation
Pn = 1.0 ./ ( alln .* alln * pi*pi/6 ) ;

% n is the number of distinct integers
n = nmax+1 ; 

% define scale of values of a for which the likelihood is required
loga = (log(1e-10):0.05:log(1e-7) ) ;
a = exp(loga) ;
sa =  size(a) ; 
numas = sa(2) ;

distbn = [0, Pn ] ;
r=1;

% Compute P(n|r) (Pngr) for all values of r and n that could be needed,
%   by explicit convolution
% Because octave's arrays start at 1, 
% row 1 is r=0; row 2 is r=1; etc.
Pngr( 1, : ) = zeros(1,n) ; 
Pngr( 1, 1 ) = 1 ;
Pngr( 2, : ) = distbn(1:n) ;
for r=2:rmax
  nextdist = conv( Pn , distbn ) ;
  distbn = [ 0,nextdist ] ;
  distbn = distbn(1:n) ;
  Pngr( r+1, : ) = distbn(1:n) ;
end
% P(n|r) made %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Make P(r|a) (Prga) for all values of r needed, and for all a
allr = (0:rmax)' ;
numr = rmax+1 ;

r = allr *  ones(1,numas) ;
aN = ones(numr,1) * a * N ;
% r and aN are both rectangular matrices
%   - rows are n; cols are a.
Prga = exp( -aN ) .* aN.**r ./ gamma(r+1) ;

% P(r|a) made %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% finally, obtain P(n|a) = sum_r P(n|r) P(r|a)
Pnga = Pngr' * Prga ;

% to get likelihood, multiply together all the required rows.
% and multiply through by 1e23 to get sensible scale
L = ones(1,numas) * 1e23 ;
sd=size(data) ;
nd = sd(2) ;
for d=1:nd
  % remember that all rows of n are offset by 1.
  L = L .* Pnga( data(d)+1 , : ) ;
end
La = [ a', L' ] ;
gset nologs y
gset logs x ;
gset nokey ;
gplot La u 1:2

ans = input ("ready")
gset size 0.4,0.4
gset yrange [1e-10:1.5] ;
gset logs y
replot