% 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