(* Package Name: salvadorV2-3.m Version: 2.3: September 10, 2008 Version: 2.2: January 28, 2008 Author: Qi Zheng, Department of Epidemiology and Biostatistics School of Rural Public Health, Texas A&M Health Science Center Starting date: March 14, 2000, National Center for Toxicological Research Major revision: March 5, 2001 Refinement started: May 11, 2001 Revisions made for Version 5.0 Mathematica: April 8, 2004 Major revision for version 2 started April 27, 2004 at School of Rural Public Health, Bryan, Texas Addition of Haldane's formulation: January 29, 2007 at SRPH Adaptation to Mathematica 6.0.1, September 18, 2007 Addition of Bartlett's formulation: December 4, 2007 Revision of Haldane's formulation added: January 20, 2008 Addition of plating efficiency: April 3, 2008 B0 limiting distribution, April 6, 2008 *) BeginPackage["salvador`"] Options[profileIntervalM]=Options[profileIntervalR]= {Alpha -> 0.05, MaxIterations -> 50, InitialM -> -1, InitialR -> -1, InitialLowerM -> -1, InitialUpperM -> -1, Tolerance -> 10.^-8, InitialLowerR -> -1, InitialUpperR -> -1, ShowIterations -> False} Options[LRIntervalLD] = {Alpha -> 0.05, Tolerance -> 10.^-8, ShowIterations->False, MaxIterations -> 50, Phi -> 1.0, InitialM->-1, InitialLowerM->-1, InitialUpperM->-1} Options[newtonLD] = {ShowIterations -> False, Phi -> 1, InitialM -> -1, MaxIterations -> 50, Tolerance -> 10.^-8} Options[newtonM] = {ShowIterations -> False, MaxIterations -> 50, Tolerance -> 10.^-8, InitialM -> -1, InitialR -> -1} Options[plotFit] = {logPlot -> False, DifferentialGrowth -> False, dataStyle -> {RGBColor[1, 0, 0], PointSize[0.015]}, fitStyle -> {RGBColor[0, 1, 0], Thickness[0.006]}, InitialM -> -1, InitialR->-1, Phi->1} (* 4-23-2005, Phi added *) Newcombe::usage="Newcombe[data] estimates m using m=(max y-mean y)/n." jonesMedian::usage="jonesMedian[data] estimates m by the method of the median of Jones et al. (1994)." fisherInfoMatrix::usage="fisherInfoMatrix[m,r,n] computes the expected Fisher information matrix for the M(m,r) model, using n terms." profileIntervalM::usage="profileIntervalM[data,opt] computes profile likelihood based confidence intervals for m under the M(m,r) model." profileIntervalR::usage="profileIntervalR[data,opt] computes profile likelihood based confidence intervals for r(=1/\!\(\[Rho]\)) under the M(m,r) model." LDLikelihood::usage="LDLikelihood[m,phi,data] gives minus log likelihood under the LD(m,phi) distribution." logLikelihoodM::usage="logLikelihoodM[m,r,data] gives the log likelihood function under the M(m,r) model." scoreInfoM::usage="scoreInfoM[data,m,r] is a temporary function ...." newtonLD::usage="newtonLD[data,opts] uses the Newton-Raphson method to compute mle of m for the LD model." newtonM::usage="newtonM[data,opts] uses the Newton-Raphson method to compute MLEs of m and r for the M(m,r) model." LRIntervalLD::usage="LRIntervalLD[X,opts] uses data X to construct a likelihood ratio based confidence interval for m that has asymptotic confidence coefficient \!\(1-\[Alpha]\); \!\(\[Alpha]\) is specified by Alpha->\!\(\[Alpha]\). Data are assumed to be generated by an LD(m,\!\(\[Phi]\)) distribution; \!\(\[Phi]\) is specified by Phi->\!\(\[Phi]\)." pmfLD::usage="pmfLD[\!\(m,\[Phi],k\)] returns a list of probabilities \!\({p\_0,p\_1,..., p\_k}\) according to a Luria-Delbr\!\(\[UDoubleDot]\)uck distribution, LD(m,\!\(\[Phi]\)). Note pmfLD[m,k] is the same as pmfLD[m,1,k]." pmfM::usage="pmfM[m,r,k] returns a list of probabilities \!\({p\_0,p\_1,...,p\_k}\) according to an \!\(M(m,r,1)\) distribution; pmfM[m,\!\(r,\[Phi]\),n] returns the same list according to an \!\(M(m,r,\[Phi])\) distribution. Note m is the mean number of mutations, \!\(r=\[Beta]\_1/\[Beta]\_2\) is the ratio of nonmutant growth rate to mutant growth rate, and \!\(\[Phi]=1-N\_0/N\_T\)." pmfMKoch::usage="pmfMKoch[m,r,k] computes \!\({p\_0,p\_1,...,p\_k}\) for an \!\(M(m,r,1)\) distribution, using Koch's algorithm. Note that r is the reciprocal of Koch's original parameter b (called \!\(\[Rho]\) in SALVADOR v.1). (This function is not recommended for large k.)" pmfPoisson::usage="pmfPoisson[\!\(\[Lambda]\),k] returns a list of probabilities \!\({p\_0,p\_1,...,p\_k}\) according to a Poisson(\!\(\[Lambda]\)) distribution." pmfCairns::usage="pmfCairns[\!\(m,\[Phi],\[Lambda],k\)] returns a list of probabilities \!\({p\_0,...,p\_k}\) according to the convolution of an LD(\!\(m,\[Phi]\)) distribution and a Poisson(\!\(\[Lambda]\)) distribution." pmf2cdf::usage="pmf2cdf[\!\({p\_0,p\_1,...,p\_k}\)] returns a list of cumulative probabilities according to the given list of probabilities \!\({p\_0,p\_1,...,p\_k}\)." simuMutants::usage="simuMutants[\!\(\[Mu],\[Beta]\_1,\[Beta]\_2, N\_0,N\_T\)] simulates the number of mutants in a sister culture. Note \!\(\[Mu]\) is the mutation rate per cell per unit time; \!\(\[Beta]\_1\) is growth rate for nonmutants; \!\(\[Beta]\_2\) is growth rate of mutants; \!\(N\_0\) is the initial population size; \!\(N\_T\) is the population size prior to plating." MeanOfMutants::usage="MeanOfMutants[\!\(\[Mu],\[Beta]\_1,\[Beta]\_2,N\_0,N\_T\)] computes the mean number of mutants on each sister culture." VarOfMutants::usage="VarOfMutants[\!\(\[Mu],\[Beta]\_1,\[Beta]\_2,N\_0,N\_T\)] computes the variance of the number of mutants on each sister culture." mleLD::usage="mleLD[data,\!\(\[Phi]:=1\),\!\(M\_\(min\)\):=0.001, \!\(M\_\(max\)\):=50] searches in the interval [\!\(M\_\(min\), M\_\(max\)\)] for the maximum likelihood estimate of m, when data is assumed to be generated by an LD(m,\!\(\[Phi]\)) distribution. Note \!\(\[Phi]\) is considered as known." mleM::usage="mleM[data,M:50,R:5] searches for the maximum likelihood estimates of both m and r by assuming that the data is generated by an M(m,r,1) distribution. The search is confined to \!\(m\[Element](0.001,M)\) and \!\(r\[Element](0.01,R)\). Note \!\(r=1/\[Rho]\)." fisherInfoLD::usage="fisherInfoLD[m,\!\(\[Phi]\),k] computes the Fisher information for m according to an LD(m,\!\(\[Phi]\)) distribution; it retains k terms in the series \!\(I(m)=\[Sum] f\^2\%j/p\_j\)." fisherInfoM::usage="fisherInfoM[m,\!\(\[Rho],\[Phi]\),k] computes the Fisher information for m, I(m), according to an M(m,\!\(\[Rho],\[Phi]\)) distribution; it retains k terms in the series \!\(I(m)=\[Sum] f\^2\%j/p\_j\). Note that fisherInfoM[m,\!\(\[Rho]\),k] assumes \!\(\[Phi]=1\)." fisherInfoRho::usage="fisherInfoRho[m,\!\(\[Rho]\),k] computes the Fisher information for \!\(\[Rho]\), I(\!\(\[Rho]\)), using an M(m, \!\(\[Rho],1\)) distribution. It retains k terms in the series \!\(I(\!\(\[Rho]\))=\[Sum] f\^2\%j/p\_j\)." fisherInfoMLea::usage="fisherInfoMLea[m,\!\(\[Rho]\),k] computes the Fisher information for both m and \!\(\[Rho]\) according to an M(m,\!\(\[Rho]\),1) distribution by the algorithm proposed by Lea and Coulson. It uses k terms to approximate I(m) and I(\!\(\[Rho]\)). (This function is not recommended for large k.)" MethodOfP0::usage="MethodOfP0[data] estimates m, the expected number of mutations, using the \!\(P\_0\) method." Options[MethodOfMeans] = {Correction -> 1, ShowIterations -> False, Tolerance -> 1.*^-8, MaxIterations -> 20, InitialM -> 1} MethodOfMeans::usage="MethodOfMeans[data,opts] iteratively solves for m the equation \!\( m log(\[Beta] n m)=X\&_ \), where \!\(\[Beta]\) is a correction factor, where \!\(X\&_\) is the sample mean, and where n is sample size. The correction suggested by Armitage corresponds to \!\(\[Beta]=3.46\)" Options[drakeMedian] = {ShowIterations -> False, InitialM->1, MaxIterations->20, Tolerance->10^-8.} (* --- July 18, 2005 --- *) drakeMedian::usage="drakeMedian[data,opts] iteratively solve the equation m log(m)=\!\(\[Xi]\&^\_0.5\), where \!\(\[Xi]\&^\_0.5\) is the sample median and where m is the expected number of mutations per culture." MethodOfMedians::usage="MethodOfMedians[data,opts] iteratively solve for m Lea and Coulson's equation \!\(\[Xi]\&^\_0.5/m\)-log(m)-1.24=0, where \!\(\[Xi]\&^\_0.5\) is the sample median. When the sample median \!\(\[Xi]\&^\_0.5\)=0, it estimate m by \!\(e\^-1.24\) which may not be a reliable estimate." Options[MethodOfMedians] = {ShowIterations -> False, InitialM->1, MaxIterations->20, Tolerance->10^-8.} semiMLE::usage="semiMLE[data,\!\(m\_0:=1\)] computes an estimate of m by solving \!\(\[Sum]Y\_i=0\), where \!\(Y\_i\) are the Lea-Coulson transforms of the original data. Note that \!\(m\_0\) is a starting value used to solve the equation." CIMutationRate::usage="CIMutationRate[data,\!\(\[Alpha],N\_T,k\)] constructs a \!\(\[Alpha]\) level confidence interval for \!\(\[Mu]\_\[Beta]\) using an LD(m,1) distribution. It uses k terms in computing the Fisher information for m. It is important to choose a suitable k. Note \!\(N\_T\) is the number of nonmutants prior to plating." plotFit::usage="plotFit[data,opts] fits either an LD(m,\!\(\[Phi]\)) or an M(m,r,1) distribution to data, and then plots the fitted distribution along with the data. Use options for different plotting styles." (* =================================================================== *) (* zero counterparts of three functions: pmfLD, newtonLD, LRIntervalLD *) (* April 21, 2006 *) Options[newtonLD0] = {ShowIterations -> False, Phi -> 1, InitialM -> -1, MaxIterations -> 50, Tolerance -> 10^-8} Options[LRIntervalLD0] = {Alpha -> 0.05, Tolerance -> 10^-8, ShowIterations->False, MaxIterations -> 50, Phi -> 1.0, InitialM->-1, InitialLowerM->-1, InitialUpperM->-1} pmfLD0::usage="pmfLD0[\!\(m,\[Phi],k\)] returns a list of probabilities \!\({p\_0,p\_1,..., p\_k}\) according to a Luria-Delbr\!\(\[UDoubleDot]\)uck distribution, LD(m,\!\(\[Phi]\)). Note pmfLD[m,k] is the same as pmfLD[m,1,k]." newtonLD0::usage="newtonLD0[data,opts] uses the Newton-Raphson method to compute mle of m for the LD model." LRIntervalLD0::usage="LRIntervalLD0[X,opts] uses data given in X to construct a likelihood ratio (LR) based confidence interval for m that has asymptotic confidence coefficient \!\(1-\[Alpha]\); \!\(\[Alpha]\) is specified by Alpha->\!\(\[Alpha]\). Data in X is assumed to be generated by an LD(m,\!\(\[Phi]\)) distribution; \!\(\[Phi]\) is specified by Phi->\!\(\[Phi]\)." (* ************* end of zero counterparts ****************** *) (* ========================================================== *) (* Haldane's formulation *) (* code development: April 27, 2006; incorporation January 29, 2007 *) (* ------------- slow -------------- *) Options[drakeHaldane] = {ShowIterations -> False, initialM->1, MaxIterations->20, Tolerance->10^-8., UseMedian->True} derivHaldane0::usage="[gen_, mu_, n_, N0_] is interesting ..." drakeHaldane::usage="drakeHaldane[data,Ng,opts] is to compute Drake's estimation." Options[newtonHaldane0] = {ShowIterations -> False, InitialMu -> -1, Tolerance -> 1.*^-8, MaxIterations -> 50, InitialCells -> 1} Options[CIHaldane0] = {Alpha -> 0.05, InitialCells -> 1, MaxIterations -> 20, Tolerance -> 1.*^-6, ShowIterations -> False, InitialLowerLimit -> -1, InitialUpperLimit -> -1, InitialMu -> -1} Options[simuHaldane] = {ShowGrowth -> False} haldaneMutations::usage="haldaneMutations[\!\(\[Mu]\),g,N0:=1] gives the expected number of mutations accumulated in the gth generation. If Kendall's Assumption C or the Galton-Waston formulation is adopted, replacing \!\(\[Mu]\) by \!\(2\[Mu]\) yields the desired result." newtonHaldane0::usage="newtonHaldane0[data,g,opts] compute the maximum likelihood estimate for \!\(\[Mu]\)." CIHaldane0::usage="CIHaldane0[data,g,opts] compute an asymptotic confidence interval for the mutation rate." rossman::usage="rossman[data,g,N0] estimates the mutation rate \!\(\[Mu]\) by the method of the mean." pmfHaldane0::usage="pmfHaldane[g, \!\(\[Mu]\), n, \!\(N\_0\)] computes the first n probabilities p(0), p(1), ..., p(n), assuming the population starts with N0 nonmutant cells. Default value for \!\(N\_0\) is unity (1)." simuHaldane::usage="simuHaldane[g, \!\(\[Mu]\)] simulates Haldane's model for g generations. The evolution of the mutation process can be viewed by setting ShowGrowth->True." (* -------------------- fast ----------------- *) Options[newtonHaldane] = {ShowIterations -> False, InitialMu -> -1, Tolerance -> 1.*^-8, MaxIterations -> 50, InitialCells -> 1} Options[CIHaldane] = {Alpha -> 0.05, InitialCells -> 1, MaxIterations -> 20, Tolerance -> 1.*^-6, ShowIterations -> False, InitialLowerLimit -> -1, InitialUpperLimit -> -1, InitialMu -> -1} Options[simuHaldane] = {ShowGrowth -> False} Options[simuAngerer] = {InitialSize -> 1, ShowGrowth -> False} haldaneMutations::usage="haldaneMutations[\!\(\[Mu]\),g,N0:=1] computes the total number of mutations accumulated in the gth generation. If Assumption C (Kdenall) or the Galton-Waston formulation is adopted, replacing \!\(\[Mu]\) by \!\(2\[Mu]\) yields the desired result." newtonHaldane::usage="newtonHaldane[data,g,opts] compute the maximum likelihood estimate for \!\(\[Mu]\)." CIHaldane::usage="CIHaldane[data,g,opts] compute an asymptotic confidence interval for the mutation rate." pmfHaldane::usage="pmfHaldane[g, \!\(\[Mu]\), n, \!\(N\_0\)] computes the first n probabilities p(0), p(1), ..., p(n), assuming the population starts with N0 nonmutant cells. Default value for \!\(N\_0\) is unity (1)." simuAngerer::usage="simuAngerer[\!\(\[Mu]\), d] simulates the first d cellular divisions using a discrete-time mutation model. Each division of a nonmutant cell can give rise to a mutant daughter cell with probability \!\(\[Mu]\). The evolution of the mutation process can be viewed by setting ShowGrowth->True, and the initial population size can be set to n0 by setting InitialSize->n0." (* ============ end of Haldane formulation ================= *) (* ========================================================== *) (* Bartlett's formulation *) (* code development: April 11, 2007; incorporation December 4, 2007 *) Options[newtonBartlett] = {InitialMutationRate -> -1, ShowIterations -> False, MaxIterations -> 20, Tolerance -> 1.*^-8} Options[CIBartlett] = {Alpha -> 0.05, InitialMutationRate -> -1, InitialLowerLimit -> -1, InitialUpperLimit -> -1, ShowIterations -> False, MaxIterations -> 20, Tolerance -> 1.*^-8} Options[goldenBartlett] = {MaxIterations -> 50, InitialTriplet -> {1, 1, 1}, Tolerance -> 0.00001, ShowIterations -> False} Options[simuBartlett] = {MaxEvents -> 1000000, ShowGrowth -> False} bartlettMean::usage="bartlettMean[X,\!\(N\_0,\[Phi]\)] estimate the mutation\ rate \!\(\[Alpha]\) by the method of moments, under the Bartlett formulation." simuBartlett::usage="simuBartlett[\!\(\[Beta],\[Mu],N\_0\),T] simulates the numbers of nonmutants and mutants using the Bartlett distribution. The mutation rate is \!\(\[Alpha]=\[Mu]/\[Beta]\) and the parameter \!\(\[Phi]\) is defined by \!\(\[Phi]=1-e\^-\[Beta]T\). The evolution of the mutation process can be viewed by setting ShowGrowth->True." pmfBli::usage="pmfBli[A,n] computes the first n+1 probabilities of a limiting distribution due to Bartlett." pmfBartlett::usage="pmfBartlett[\!\(\[Alpha],\[Phi]\),n] computes the first n+1 probabilities of a Bartlett distribution having \!\(N\_0=1\). Note that pmfBartlett[\!\(\[Alpha],\[Phi]\),n,k] assumes \!\(N\_0\)=k." newtonBartlett::usage="newtonBartlett[X, \!\(N\_0, \[Phi]\)] finds the maximum\ likelihood estimate of the mutation rate \!\(\[Alpha]\) under the Bartlett\ formulation." CIBartlett::usage="CIBartlett[X, \!\(N\_0, \[Phi]\)] constructs an asymptotic confidence interval for the mutation rate \!\(\[Alpha]\) under the Bartlett formulation." goldenBartlett::usage="goldBartlett[X, \!\(N\_0, \[Phi]\)] finds the maximum\ likelihood estimate of the mutation rate by the golden section search method,\ under the Bartlett formulation." (* --------- zero functions ---------- *) Options[newtonBartlett0] = {InitialMutationRate -> -1, ShowIterations -> False, MaxIterations -> 20, Tolerance -> 1/100000000} Options[CIBartlett0] = {Alpha -> 5/100, InitialMutationRate -> -1, InitialLowerLimit -> -1, InitialUpperLimit -> -1, ShowIterations -> False, MaxIterations -> 20, Tolerance -> 1/100000000} pmfBartlett0::usage="pmfBartlett[\!\(\[Alpha],\[Phi]\),n] computes the first\ n probabilities of a Bartlett distribution. Note that\ pmfBartlett[\!\(\[Alpha],\[Phi]\),n,k] will assume \!\(N\_0\)=k." newtonBartlett0::usage="newtonBartlett[X, \!\(N\_0, \[Phi]\)] finds the maximum\ likelihood estimate of the mutation rate \!\(\[Alpha]\) under the Bartlett\ formulation." CIBartlett0::usage="CIBartlett0[X, \!\(N\_0, \[Phi]\)] finds a confidence\ interval for the mutation rate \!\(\[Alpha]\) under the Bartlett formulation." listConv::usage="listConv[x,y] does convolution of x and y." goldenBartlett0::usage="goldBartlett0[X, \!\(N\_0, \[Phi]\)] finds the maximum\ likelihood estimate of the mutation rate by the golden section search method." (* ============ end of Bartlett's formulation ================= *) (* ============ Plating efficiency, April 3, 2008 =============== *) Options[newtonLDPlating] = {InitialM -> -1, ShowIterations -> False, MaxIterations -> 20, Tolerance -> 1.0*^-8} Options[empiricalMLE] = {InitialM -> 1, ShowIterations -> False, MaxIterations -> 20, Tolerance -> 1.*^-8} Options[CILDPlating] = {Alpha -> 0.05, InitialM -> -1, InitialLowerLimit -> -1, InitialUpperLimit -> -1, ShowIterations -> False, MaxIterations -> 20, Tolerance -> 1.*^-8} Options[LCAMedian] = {InitialM -> 1, ShowIterations -> False, MaxIterations -> 20, Tolerance -> 1.*^-8} etaSeq::usage="etaSeq[e,n] yields the eta sequence ..., March 12, 2008." LDLoglikely$Score4Plating::usage="LDLoglikely$Score4Plating[m,e,data] is temporary ..." empiricalMLE::usage="empricalMLE[data,e,opts] estimates m when plating efficiency is e." Nadas::usage="Nadas[data,Nt,N0] computes the mutation rate by assuming synchronized growth." LCAMedian::usage="LCAMedian[data,e,opts] estimates m according to a formula of Lea and Coulson modified by Angerer." mleAngerer::usage="mleAngerer[data,e,guess:1] computes MLE of m by a psudo maximum likelihood approach, where e is known plating efficiency." newtonLDPlating::usage="newtonLDPlating[data,\!\(\[Epsilon]\),opts] estimates the expected number of mutations per culture under the assumption of imperfect plating efficiency, that is, \!\(0<\[Epsilon]<1\)." newtonLDPlating0::usage="newtonLDPlating0[data,\!\(\[Epsilon]\),opts] estimates the expected number of mutations per culture under the assumption of imperfect plating efficiency, that is, \!\(0<\[Epsilon]<1\)." CILDPlating::usage="CILDPlating[data,\!\(\[Epsilon]\),opts] constructs a \!\((1-\[Alpha])\[RawPercent]\) confidence interval for m, the expected number of mutations per culture. Plating efficiency \!\(\[Epsilon]\) must satisfy \!\(0<\[Epsilon]<1\). " CILDPlating0::usage="CILDPlating0[data,\!\(\[Epsilon]\),opts] constructs a \!\((1-\[Alpha])\[RawPercent]\) confidence interval for m, the expected number of mutations per culture. Plating efficiency \!\(\[Epsilon]\) must satisfy \!\(0<\[Epsilon]<1\). " P0Plating::usage="P0Plating[data,\!\(\[Epsilon]\)] estimate the expected number of mutations per culture using the P0 method. It is assumed that the plating efficiency \!\(\[Epsilon]\) satisfies 0<\!\(\[Epsilon]\)<1." JonesMedianPlating::usage="JonesMedianPlating[data,\!\(\[Epsilon]\)] estimate the expected number of mutations per culture using the method of the median of Jones et al. It is assumed that the plating efficiency \!\(\[Epsilon]\) satisfies 0<\!\(\[Epsilon]\)<1." pmfLDPlating::usage="pmfLDPlating[m,\!\(\[Epsilon]\),n] computes probabilities for an LD(m) distribution with plating efficiency \!\(\[Epsilon]\)." pmfLDPlating0::usage="pmfLDPlatingF[m,\!\(\[Epsilon]\),n] computes probabilities for an LD(m) distribution with plating efficiency \!\(\[Epsilon]\)." (* ========= end of plating ========= *) (* ************ B0 distribution ************ *) Options[newtonBartzero] = {InitialA -> -1, ShowIterations -> False, MaxIterations -> 20, Tolerance -> 1.*^-6} Options[CIBartzero] = {Alpha -> 0.05, InitialA -> -1, InitialLowerLimit -> -1, InitialUpperLimit -> -1, ShowIterations -> False, MaxIterations -> 20, Tolerance -> 1.*^-8} p0Bartzero::usage="p0Bartzero[data_, k_] estimates the parameter A using the P0 methods when k is known." simuBartzero::usage="simuBartzero[A,k,n,M] simulate n copies of the \!\(B\^0(A,k)\) distribution, using a look-up table of length M." newtonBartzero::usage="newtonBartzero[data,k,opts] used Newton-Raphson ..." CIBartzero::usage="CIBartzero[data,k,opts] assumes k is know, and finds mle of A." pmfBartzero::usage="pmfBartzero[A,k, n] computes the first n+1 probabilities of the Bartlett limiting distribution. Note the parameter k can be a fractional." (* ====== end of B0 distribution ======== *) Begin["`Private`"] (* Fisher's information I(m) for an LD(m,phi) *) (* revised October 30, 2001 *) fisherInfoLD[m_, phi_:1, n_] := Module[{l, p, f}, l = Table[phi^(k-1)*(1/k - phi/(k + 1)), {k, 1, n}]; l = Prepend[l, -1]; p = pmfLD[m, phi, n]; If[Last[p] == 0, Return[Print["n is too large \!\(p\_n=0\) encountered."]]]; f = $luria`$pdfconv[p, N[l]]; Plus @@ (f^2/p)]/;(m>0)&&(phi>0)&&(phi<=1)&&(n>0) (* observed information J(m) for an LD(m,phi) 4/27/04 *) ObservedInfoLD[m_, phi_:1, U_] := Module[{l, p, p1, p2}, l = Table[phi^(k-1)*(1/k - phi/(k + 1)), {k, 1, U}]; l = N[Prepend[l, -1]]; p = pmfLD[m, phi, U]; p1=$luria`$pdfconv[p,l]; p2=$luria`$pdfconv[p1,l]; Plus@@(p1^2/p-p2) ]/;(m>0)&&(phi>0)&&(phi<=1)&&(U>0) crossInfo[m_, r_, n_] := Module[{Lm, Lr, be, et, p, p10, p01, f11, f12, f22}, be = Table[(j - 1)/(j + r), {j, 2, n}]; be = FoldList[Times, 1, be]; be = be/(1 + r); Lm = Prepend[r*be, -1]; et = Table[1/(j + r), {j, 1, n}]; et = FoldList[Plus, First[et], Rest[et]]; Lr = m*be*(1 - r*et); Lr = Prepend[Lr, 0]; p = pmfM[m, r, 1, n]; p10=$luria`$pdfconv[p,Lm]; p01=$luria`$pdfconv[p,Lr]; f11=p10*p10/p; f12=p10*p01/p; f22=p01*p01/p; {Plus@@f11, Plus@@f12, Plus@@f22}] (* Fisher's information I(m) for an M(m,rho,phi) distribution *) fisherInfoM[m, 1, phi_, n_Integer]:=fisherInfoLD[m,phi,n] fisherInfoM[m_, r_, 1, n_Integer] := $luria`$kochinfo[N[m],N[r],n] fisherInfoM[m_, r_, n_Integer] := fisherInfoM[m,r,1,n] fisherInfoM[m_, r_, phi_, n_Integer] := Module[{l, p, f, omega}, omega = 1 - (1 - phi)^r; l = (Beta[omega, #1, 1 + 1/r] & ) /@ Range[n]; l = Prepend[l/r/phi, -1]; p = pmfM[m, r, phi, n]; f = $luria`$pdfconv[l, p]; Plus @@ (f^2/p)]/;(m>0)&&(r>0)&&(phi>0)&&(phi<1)&&(n>0) fisherInfoRho[m_, r_, n_Integer] := $luria`$iofrho[N[m],N[r],n] (* FisherInfoLDLC[m_/;Positive[m],n_Integer]/;n>0:= $luria`$fishinfo[N[m],n]*) fisherInfoMLea[m_,r_,n_]:=$luria`$fisherho[N[m],N[r],n] (* simulating the number of mutants in a sister culture *) simuMutants[mu_, b1_, b2_, N0_, Nt_] := Module[{meanMutation, noOfMutation, T, mutationEpoch}, T = Log[Nt/N0]/b1; meanMutation = mu*N0*(Exp[b1*T] - 1)/b1; noOfMutation = Random[PoissonDistribution[meanMutation]]; mutationEpoch = Table[Log[1 + Random[]*(Exp[b1*T] - 1)]/b1,{noOfMutation}]; Plus @@ (Random[GeometricDistribution[Exp[-b2*(T - #1)]]] + 1 & ) /@ mutationEpoch] (* mean number of mutants *) MeanOfMutants[mu_, b1_, b2_, N0_, Nt_] := Module[{t}, t = Log[Nt/N0]/b1; If[b1 == b2, mu*N0*t*Exp[b1*t], mu*N0*(Exp[b1*t] - Exp[b2*t])/(b1 - b2)]] (* variance of the number of mutants *) VarOfMutants[mu_, b1_, b2_, N0_, Nt_] := Module[{t}, t = Log[Nt/N0]/b1; If[b1 == b2, 2*mu*N0*Exp[b1*t]*(Exp[b1*t] - 1)/b1 - mu*N0*t*Exp[b1*t], If[b1 == 2*b2, mu*N0*Exp[b2*t]*(1 - Exp[b2*t] + 2*b2*t*Exp[b2*t])/b2, mu*N0*Exp[b2*t]*(b1*(1 + Exp[(b1 - b2)*t] - 2*Exp[b2*t]) + 2*b2*(Exp[b2*t] - 1))/(b1 - b2)/(b1 - 2*b2)]]] (* Poisson distribution *) pmfPoisson[lambda_ /; Positive[lambda], k_Integer /; NonNegative[k]] := Table[PDF[PoissonDistribution[lambda],j],{j,0,k}] (* Luria-Delbruck distribution *) pmfLD[m_/;Positive[m], phi_/;(phi>0)&&(phi<=1), k_Integer/;NonNegative[k]] := $luria`$pdfluria[N[m],N[phi],k] pmfLD[m_, k_Integer]:=pmfLD[m,1,k]/;(m>0)&&(k>=0) (* Cairnsian distribution *) pmfCairns[m_, phi_, poi_, k_Integer] := Module[{LDProb, PoiProb}, LDProb = $luria`$pdfluria[N[m], N[phi], k]; PoiProb = pmfPoisson[N[poi],k]; $luria`$pdfconv[LDProb,PoiProb] ]/;(m>0)&&(phi>0)&&(phi<=1)&&(poi>0)&&(k>=0) (* the M distribution *) pmfM[m_,r__,0]:={Exp[-N[m]]} (* Koch's algorithm for the M(m,r,1) distribution, slow for large n *) pmfMKoch[m_, r_, 0]:={Exp[-m]} pmfMKoch[m_, r_, n_Integer] := Module[{q}, q=$luria`$pdfkoch[N[m],N[1/r],n]; (* Koch uses rho=1/r, 5/14/04 *) Exp[-N[m]] Prepend[q,1] ]/;(m>0)&&(r>0)&&(n>0) (* pmf for an M(m,r,1) distribution, using an algorithm like in Stewart et al. *) pmfM[m_,r_,n_Integer]:=$luria`$kochstew[N[m],N[r],n] (* from 1/r to r 5/13/04 *) pmfM[m_,r_,1,n_Integer]:=$luria`$kochstew[N[m],N[r],n] (* from 1/r to r, 5/13/04 *) (* exact algorithm of Stewart et al. *) pmfM[m_, r_, phi_, n_Integer] := Module[ {q, r1=1+r, omega=N[1-(1-phi)^(1/r)]}, q = (Beta[omega, #1, r1] & ) /@ Range[n]; q= q*m*r/phi; q = Prepend[q, -N[m]]; $luria`$poisum[q] ]/;(phi>0)&&(phi<1)&&(m>0)&&(r>0)&&(n>0) (* cumulative probabilities *) pmf2cdf[l_List] := FoldList[Plus, First[l], Rest[l]] /; Length[l] > 0 (* the median method of Jones et al. 5/19/04 *) jonesMedian[data_]:=Module[{r}, r=Median[data]; (r-0.693147)/(Log[r]+0.3665)] (* the Newcombe method, 5/20/04 *) Newcombe[data_]:=Module[{},(Max[data]-Mean[data])/Length[data]//N] (* the P0 method *) MethodOfP0[data_] := Module[{p0}, p0 = Count[data, 0]/Length[data]; If[p0 > 0, -Log[N[p0]], Return[ Print[ "The \!\(P\_0 \) method is not applicable to this case."];Null]] ] (* modified method of the mean, with Armitage correction, July 26, 2005 *) MethodOfMeans[data_, opts___] := Module[ {i, med, m0, m1, showIter, guess, tol, maxIter, cor}, showIter = ShowIterations /. {opts} /. Options[MethodOfMeans]; guess = InitialM /. {opts} /. Options[MethodOfMeans]; tol = Tolerance /. {opts} /. Options[MethodOfMeans]; maxIter = MaxIterations /. {opts} /. Options[MethodOfMeans]; cor = Correction /. {opts} /. Options[MethodOfMeans]; If[guess <= 0, m0 = 1., m0 = N[guess]]; If[showIter, Print[m0]]; For[i = 1, i <= maxIter, i++, m1 = m0 - (m0*Log[m0] - cor*Length[data]*Mean[data])/(1 + Log[m0]); If[showIter, Print[m1/cor/Length[data]]]; If[Abs[(m0 - m1)/m0] < tol, Return[m1/cor/Length[data]], m0 = m1]]; Print["no convergence"]; ] (* ----------- Drake's Median Estimator ----------- July 27, 2005 *) drakeMedian[data_, opts___] := Module[ {i, med, m0, m1, showIter, guess, tol, maxIter}, showIter = ShowIterations /. {opts} /. Options[drakeMedian]; guess = InitialM /. {opts} /. Options[drakeMedian]; tol = Tolerance /. {opts} /. Options[drakeMedian]; maxIter = MaxIterations /. {opts} /. Options[drakeMedian]; med = Median[data]; If[med==0, Print["Sample median is zero, Drake's method not applicable."]; Return[]]; If[guess<=0, m0=1.0, m0 = N[guess]]; If[showIter, Print[m0]]; For[i = 1, i <= maxIter, i++, m1 = m0 - (m0*Log[m0]-med)/(1+Log[m0]); If[showIter, Print[m1]]; If[Abs[(m0 - m1)/m0] < tol, Return[m1], m0 = m1]]; Print["no convergence"];] (* --------- to recode the method of the median, January 31, 2007 ---------- *) MethodOfMedians[data_?NumericQ, opts___] := MethodOfMedians[{data},opts] MethodOfMedians[data_List, opts___] := Module[ {i, med, c, w0, w1, showIter, guess, tol, maxIter}, showIter = ShowIterations /. {opts} /. Options[MethodOfMedians]; guess = InitialM /. {opts} /. Options[MethodOfMedians]; tol = Tolerance /. {opts} /. Options[MethodOfMedians]; maxIter = MaxIterations /. {opts} /. Options[MethodOfMedians]; med = Median[data]; c=Log[med]+1.24; If[med==0, Return[Exp[-1.24]]]; If[guess<=0, w0=1.0, w0 = med/N[guess]]; If[showIter, Print[med/w0]]; For[i = 1, i <= maxIter, i++, w1 = w0 - (w0+Log[w0]-c)/(1+1/w0); If[showIter, Print[med/w1]]; If[Abs[(w0 - w1)/w0] < tol, Return[med/w1], w0 = w1]]; Print["no convergence"];] (* the semi-MLE of Lea and Coulson *) semiMLE[data_,m0_:1] := Module[{m,root}, root = FindRoot[leaLike[m, data], {m, m0}]; root[[1]][[2]] ] (* MLE for m using an LD(m,phi) distribution, phi is assumed to be known *) mleLD[data_, phi_:1, minM_:0.001, maxM_:50] := Module[{m, opt}, opt = FindMinimum[LDLikelihood[m, phi, data], {m, minM, maxM}]; opt[[2]][[1]][[2]] ] (* MLE for both m and r, using an M(m,r,1) distribution *) mleM[data_, M_:50, R_:5] := Module[{m, r, opt}, opt = FindMinimum[-logLikelihoodM[m, r, data], {{m, 0.001, M}, {r, 0.01, R}}]; {opt[[2]][[1]][[2]], opt[[2]][[2]][[2]]}] (* confidence interval for mutation rate with an LD(m,1) model *) CIMutationRate[data_, level_, Nt_, k_:1000] := Module[{m, z, d, f0, f1}, z = Quantile[NormalDistribution[0, 1], (1 - level)/2]; m = mleLD[data]; f0 = fisherInfoLD[m, 1, k]; f1 = fisherInfoLD[m, 1, k + 1000]; If[f1 - f0 > 0.001, Return[Print["Specified k is too small."]]]; d = z/Sqrt[Length[data]*f1]; {m + d, m - d}/Nt] (* plot the fitting, May 11, 2001 *) plotFit[data_List, opts___] := Module[{distinctPt, freq, cumuFreq, g1, g2, m, r, prob, koch, pointStyle, lineStyle, useLog, rescale, initM, initR, phi}, koch = DifferentialGrowth /. {opts} /. Options[plotFit]; pointStyle = dataStyle /. {opts} /. Options[plotFit]; lineStyle = fitStyle /. {opts} /. Options[plotFit]; useLog = logPlot /. {opts} /. Options[plotFit]; initM = InitialM /. {opts} /. Options[plotFit]; initR = InitialR /. {opts} /. Options[plotFit]; (* 5-17-2004 *) phi = Phi /. {opts} /. Options[plotFit]; (* April 23, 2005 *) rescale = If[useLog, (Log[#+1])&, Identity]; (* add "1" ... 5/18/04 *) distinctPt = Sort[Union[data]]; freq = (Count[data, #1] & ) /@ distinctPt; cumuFreq = FoldList[Plus, First[freq], Rest[freq]]/Length[data]; g1 = ListPlot[Transpose[{rescale[distinctPt], cumuFreq}], PlotStyle -> Flatten[{pointStyle}], DisplayFunction -> Identity]; If[koch, {m, r} = newtonM[data, InitialM->initM, InitialR->initR], m = newtonLD[data, InitialM->initM, Phi->phi]]; (* 4-30-04, Apr-23-05 *) If[koch, prob = pmfM[m, r, Max[data]], prob = pmfLD[m, 1, Max[data]]]; prob = (pmf2cdf[prob])[[Range[1+Min[data], 1+Max[data]]]]; g2 = ListPlot[Transpose[{rescale[Range[Min[data], Max[data]]], prob}], Joined -> True, PlotStyle -> Flatten[{lineStyle}], PlotRange->All, (* this option added 9-20-07 for Mathematica 6.0 *) DisplayFunction -> Identity]; Show[g1, g2, PlotRange -> All, Frame -> True, DisplayFunction -> $DisplayFunction, Axes -> None]] (* likelihood ratio-based confidence interval for m with LD distribution, 5/4/04 *) LRIntervalLD[data_List, opts___] := Module[ {mhat, la, qa, h, m0, m1, like, score, info, mL, mU, alpha, tol, initM, initLow, initUp, maxIter, phi, showIter, i}, alpha = Alpha /. {opts} /. Options[LRIntervalLD]; tol = Tolerance /. {opts} /. Options[LRIntervalLD]; maxIter = MaxIterations /. {opts} /. Options[LRIntervalLD]; phi = Phi /. {opts} /. Options[LRIntervalLD]; initM= InitialM /. {opts} /. Options[LRIntervalLD]; initLow= InitialLowerM /. {opts} /. Options[LRIntervalLD]; initUp= InitialUpperM /. {opts} /. Options[LRIntervalLD]; showIter= ShowIterations /. {opts} /. Options[LRIntervalLD]; If[phi<0||phi>1, Return["\!\(\[Phi]\) out of range!"]]; mhat = newtonLD[data, Phi->phi, InitialM->initM, Tolerance->tol]; {like, score} = LDLoglikeScore[mhat, phi, data]; {score, info} =LDScore$Info[mhat,phi,data]; qa=Quantile[ChiSquareDistribution[1], 1-alpha]; h=Sqrt[qa/info]; la = like - 0.5*qa; If[showIter, Print["Iterating for lower limit ..."]]; If[initLow<=0, m0=mhat-0.5*h, m0=initLow]; For[i = 1, i <= maxIter, i++, {like, score} = LDLoglikeScore[m0, phi, data]; m1 = m0 - (like - la)/score; If[showIter, Print[m0]]; If[Abs[(m1 - m0)/m0] <= tol, mL=m1; Break[], m0 = m1]]; If[showIter, Print["Iterating for upper limit ..."]]; If[initUp<=0, m0=mhat+0.5*h, m0=initUp]; For[i = 1, i <= maxIter, i++, {like, score} = LDLoglikeScore[m0, phi, data]; m1 = m0 - (like - la)/score; If[showIter, Print[m0]]; If[Abs[(m1 - m0)/m0] <= tol, mU=m1; Break[], m0 = m1]]; {mL, mU}] (* computing MLE of m using Newton-Raphson method, for LD model, 5/5/04 *) newtonLD[data_List, opts___] := Module[ {m0, m1, score, info, showIter, phi, tol, maxIter, initEst, i}, showIter = ShowIterations /. {opts} /. Options[newtonLD]; initEst = InitialM /. {opts} /. Options[newtonLD]; maxIter = MaxIterations /. {opts} /. Options[newtonLD]; phi = Phi /. {opts} /. Options[newtonLD]; tol = Tolerance /. {opts} /. Options[newtonLD]; If[initEst>0, m0=initEst, If[Median[data]==0, m0=MethodOfP0[data], m0=jonesMedian[data]]]; If[showIter, Print[m0]]; For[i = 1, i <= maxIter, i++, {score, info} = LDScore$Info[m0, phi, data]; m1 = m0 + score/info; If[Abs[(m1 - m0)/m0] < tol, Return[m1], If[showIter, Print[m1]]; m0 = m1]]; Print["no convergence!"]; ] (* Newton-Raphson for MLE of the M(m,r) model, 5/13/2004 *) newtonM[data_, opts___] := Module[ {showIter, maxIter, tol, init4m, init4r, m0, r0, m1, r1, score, info, i}, showIter = ShowIterations /. {opts} /. Options[newtonM]; maxIter = MaxIterations /. {opts} /. Options[newtonM]; tol = Tolerance /. {opts} /. Options[newtonM]; init4m = InitialM /. {opts} /. Options[newtonM]; init4r = InitialR /. {opts} /. Options[newtonM]; If[init4m < 0, m0 = newtonLD[data], m0 = init4m]; If[init4r < 0, r0 = 1, r0 = init4r]; {score, info} = scoreInfoM[data, N[m0], N[r0]]; For[i = 1, i <= maxIter, i++, {m1, r1} = {m0, r0} + score . Inverse[info]; If[Abs[(m1 - m0)/m0] + Abs[(r1 - r0)/r0] < tol, Return[{m1, r1}], If[showIter, Print[{m0, r0}]]; m0 = m1; r0 = r1; {score, info} = scoreInfoM[data, m1, r1]; ]; ]; Print["no convergence"]] (* confidence interval for m based on profile likelihood, for M(m,r) model *) (* 5/15/2004 *) profileIntervalM[data_, opts___] := Module[ {showIter, initM, initR, initLowM, initUpM, initLowR, initUpR, tol, alpha, maxIter, m0, m1, r0, r1, la, mhat, rhat, score, info, like, lowCL, upCL, i, venLowM, venLowR, venUpM, venUpR}, alpha = Alpha /. {opts} /. Options[profileIntervalM]; showIter = ShowIterations /. {opts} /. Options[profileIntervalM]; maxIter = MaxIterations /. {opts} /. Options[profileIntervalM]; tol = Tolerance /. {opts} /. Options[profileIntervalM]; initM = InitialM /. {opts} /. Options[profileIntervalM]; initLowM = InitialLowerM /. {opts} /. Options[profileIntervalM]; initUpM = InitialUpperM /. {opts} /. Options[profileIntervalM]; initR = InitialR /. {opts} /. Options[profileIntervalM]; initLowR = InitialLowerR /. {opts} /. Options[profileIntervalM]; initUpR = InitialUpperR /. {opts} /. Options[profileIntervalM]; {mhat, rhat} = newtonM[data, InitialM -> initM, InitialR -> initR, Tolerance -> tol]; la = logLikelihoodM[mhat, rhat, data] - 0.5*Quantile[ChiSquareDistribution[1], 1 - alpha]; {venLowM, venLowR, venUpM, venUpR}=venzon4m[data,mhat,rhat,alpha]; If[initLowM < 0, m0 = venLowM, m0 = initLowM]; If[initLowR < 0, r0 = venLowR, r0 = initLowR]; If[showIter, Print["Iterating for lower limit ..."]]; {score, info} = scoreInfoM[data, m0, r0]; like = logLikelihoodM[m0, r0, data]; For[i = 1, i <= maxIter, i++, {m1, r1} = {m0, r0} + Inverse[{-score, {info[[1,2]], info[[2,2]]}}] . {like - la, score[[2]]}; If[showIter, Print[{m0, r0}]]; {score, info} = scoreInfoM[data, m1, r1]; If[Abs[(m0 - m1)/m0] < tol, lowCL = m1; Break[]]; like = logLikelihoodM[m1, r1, data]; m0 = m1; r0 = r1; ]; (* upper confidence limit *) If[initUpM < 0, m0 = venUpM, m0 = initUpM]; If[initUpR < 0, r0 = venUpR, r0 = initUpR]; If[showIter, Print["Iterating for upper limit ..."]]; {score, info} = scoreInfoM[data, m0, r0]; like = logLikelihoodM[m0, r0, data]; For[i = 1, i <= maxIter, i++, {m1, r1} = {m0, r0} + Inverse[{-score, {info[[1,2]], info[[2,2]]}}] . {like - la, score[[2]]}; If[showIter, Print[{m0, r0}]]; {score, info} = scoreInfoM[data, m1, r1]; If[Abs[(m0 - m1)/m0] < tol, upCL = m1; Break[]]; like = logLikelihoodM[m1, r1, data]; m0 = m1; r0 = r1; ]; {lowCL, upCL}] (* profile interval for r=1/rho under the M(m,r) model, 5/18/04 *) profileIntervalR[data_, opts___] := Module[ {showIter, initM, initR, initLowM, initUpM, initLowR, initUpR, tol, alpha, maxIter, m0, m1, r0, r1, la, mhat, rhat, score, info, like, lowCL, upCL, i, venLowM, venLowR, venUpM, venUpR}, alpha = Alpha /. {opts} /. Options[profileIntervalR]; showIter = ShowIterations /. {opts} /. Options[profileIntervalR]; maxIter = MaxIterations /. {opts} /. Options[profileIntervalR]; tol = Tolerance /. {opts} /. Options[profileIntervalR]; initM = InitialM /. {opts} /. Options[profileIntervalR]; initLowM = InitialLowerM /. {opts} /. Options[profileIntervalR]; initUpM = InitialUpperM /. {opts} /. Options[profileIntervalR]; initR = InitialR /. {opts} /. Options[profileIntervalR]; initLowR = InitialLowerR /. {opts} /. Options[profileIntervalR]; initUpR = InitialUpperR /. {opts} /. Options[profileIntervalR]; {mhat, rhat} = newtonM[data, InitialM -> initM, InitialR -> initR, Tolerance -> tol]; la = logLikelihoodM[mhat, rhat, data] - 0.5*Quantile[ChiSquareDistribution[1], 1 - alpha]; {venLowM, venLowR, venUpM, venUpR}=venzon4r[data,mhat,rhat,alpha]; If[initLowM < 0, m0 = venLowM, m0 = initLowM]; If[initLowR < 0, r0 = venLowR, r0 = initLowR]; If[showIter, Print["Iterating for lower limit ..."]]; {score, info} = scoreInfoM[data, m0, r0]; like = logLikelihoodM[m0, r0, data]; For[i = 1, i <= maxIter, i++, {m1, r1} = {m0, r0} + Inverse[{-score, {info[[1,1]], info[[1,2]]}}] . {like - la, score[[1]]}; If[showIter, Print[{m0, r0}]]; {score, info} = scoreInfoM[data, m1, r1]; If[Abs[(r0 - r1)/r0] < tol, lowCL = r1; Break[]]; like = logLikelihoodM[m1, r1, data]; m0 = m1; r0 = r1; ]; (* upper confidence limit *) If[initUpM < 0, m0 = venUpM, m0 = initUpM]; If[initUpR < 0, r0 = venUpR, r0 = initUpR]; If[showIter, Print["Iterating for upper limit ..."]]; {score, info} = scoreInfoM[data, m0, r0]; like = logLikelihoodM[m0, r0, data]; For[i = 1, i <= maxIter, i++, {m1, r1} = {m0, r0} + Inverse[{-score, {info[[1,1]], info[[1,2]]}}] . {like - la, score[[1]]}; If[showIter, Print[{m0, r0}]]; {score, info} = scoreInfoM[data, m1, r1]; If[Abs[(r0 - r1)/m0] < tol, upCL = r1; Break[]]; like = logLikelihoodM[m1, r1, data]; m0 = m1; r0 = r1; ]; {lowCL, upCL}] (* expected Fisher Information matrix for the M(m,r) model, 5/17/2004 *) fisherInfoMatrix[m_?NumericQ, r_?NumericQ, n_] := Module[ {B, eta, h10, h01, h02, p, p10, p01, I11, I12, I22}, B = N[Table[(j - 1)/(j + r), {j, 2, n}]]; B = FoldList[Times, 1, B]/(1 + r); eta = Table[1/(j + r), {j, 1, n}]; eta = FoldList[Plus, First[eta], Rest[eta]]; h10 = Prepend[r*B, -1]; h01 = Prepend[B*(1 - r*eta), 0]; p = $luria`$kochstew[N[m], N[r], n]; p10 = $luria`$pdfconv[p, h10]; p01 = m*$luria`$pdfconv[p, h01]; I11 = Plus@@ (p10^2/p); I12 = Plus@@ (p10*p01/p); I22 = Plus@@ (p01^2/p); {{I11, I12}, {I12, I22}}] (* ============== internal functions ============== *) LDLikelihood[m_?NumericQ,phi_,data_List] := Module[{prob, likelihood}, prob = $luria`$pdfluria[N[m], N[phi], Max[data]]; likelihood = (-Log[prob[[#1 + 1]]] & ) /@ data; Plus @@ likelihood] (* loglikelihood and score, 5-3-04 *) LDLoglikeScore[m_?NumericQ,phi_,data_List] := Module[ {n, h, p, p1, loglike, score}, n= Max[data]; p = $luria`$pdfluria[N[m], N[phi], n]; loglike = Log[p[[#1 + 1]]] & /@ data; h = Table[phi^(k-1)*(1/k - phi/(k + 1)), {k, 1, n}]; h = Prepend[h, -1]; p1 = $luria`$pdfconv[p, N[h]]; score = ((p1/p)[[#1 + 1]])& /@ data; {Plus @@ loglike, Plus @@ score} ] (* Rao's score U and observed Fisher information J for LD model, 5-5-04 *) LDScore$Info[m_?NumericQ, phi_, data_List] := Module[ {n, h, p, p1, p2, info, score}, n = Max[data]; p = $luria`$pdfluria[N[m], N[phi], n]; h = Table[phi^(k - 1)*(1/k - phi/(k + 1)), {k, 1, n}]; h = Prepend[h, -1]; p1 = $luria`$pdfconv[p, N[h]]; p2 = $luria`$pdfconv[p1, N[h]]; score = ((p1/p)[[#1 + 1]] & ) /@ data; info = (((p1/p)^2 - p2/p)[[#1 + 1]] & ) /@ data; {Plus @@ score, Plus @@ info}] (* all 1st and 2nd derivatives of the log likelihood function for M(m,r) model *) (* 5/13/2004 *) partials4M[m_, r_, n_] := Module[ {B, eta, h10, h01, h02, p, p10, p01, p02, p20, p11}, B = Table[(j - 1)/(j + r), {j, 2, n}]; B = FoldList[Times, 1, B]/(1 + r); eta = Table[1/(j + r), {j, 1, n}]; eta = FoldList[Plus, First[eta], Rest[eta]]; h10 = Prepend[r*B, -1]; h01 = Prepend[B*(1 - r*eta), 0]; p = $luria`$kochstew[N[m], N[r], n]; p10 = $luria`$pdfconv[p, h10]; p01 = m*$luria`$pdfconv[p, h01]; h02 = Table[j/(j + r)^2, {j, 1, n}]; h02 = FoldList[Plus, First[h02], Rest[h02]]; h02 = B*(r*eta*eta - eta - h02); h02 = Prepend[h02, 0]; p02 = m*($luria`$pdfconv[h02, p] + $luria`$pdfconv[h01, p01]); p20 = $luria`$pdfconv[h10, p10]; p11 = (1*p01)/m + m*$luria`$pdfconv[h01, p10]; {p10, p01, p20, p11, p02, p}] (* score and observed information matrix for the M(m,r) model, 5/13/04 *) scoreInfoM[data_List, m_, r_] := Module[ {U1, U2, J11, J12, J22, p10, p01, p20, p11, p02, p, n}, n = Max[data]; {p10, p01, p20, p11, p02, p} = partials4M[m, r, n]; U1 = ((p10/p)[[#1 + 1]] & ) /@ data; U1 = Plus @@ U1; U2 = ((p01/p)[[#1 + 1]] & ) /@ data; U2 = Plus @@ U2; J11 = (((p10/p)^2 - p20/p)[[#1 + 1]] & ) /@ data; J11 = Plus @@ J11; J12 = (((p10*p01)/p^2 - p11/p)[[#1 + 1]] & ) /@ data; J12 = Plus @@ J12; J22 = (((p01/p)^2 - p02/p)[[#1 + 1]] & ) /@ data; J22 = Plus @@ J22; {{U1, U2}, {{J11, J12}, {J12, J22}}}] logLikelihoodM[m_?NumericQ,r_,data_List] := Module[{prob, likelihood}, prob = pmfM[N[m], N[r], Max[data]]; likelihood = (Log[prob[[#1 + 1]]] & ) /@ data; Plus @@ likelihood] (* adapted from MLikelihood, 5/14/04 *) (* find initial estimate of confidence limits for m and r according to Venzon *) (* output give a list of the form {lower m, lower r, upper m, upper r} *) (* see Venzon and Moolgavkar 1988, 5/15/04 *) venzon4m[data_, mhat_, rhat_, alpha_:0.05] := Module[{U, J, q, h}, {U, J} = scoreInfoM[data, mhat, rhat]; q = Quantile[ChiSquareDistribution[1], 1 - alpha]; h = J[[1,1]] - J[[1,2]]^2/J[[2,2]]; h = 0.5*Sqrt[q/h]; {mhat - h, rhat + (h*J[[1,2]])/J[[2,2]], mhat + h, rhat - (h*J[[1,2]])/J[[2,2]]} ] venzon4r[data_, mhat_, rhat_, alpha_:0.05] := Module[{U, J, q, h}, {U, J} = scoreInfoM[data, mhat, rhat]; q = Quantile[ChiSquareDistribution[1], 1 - alpha]; h = J[[2,2]] - J[[1,2]]^2/J[[1,1]]; h = 0.5*Sqrt[q/h]; {mhat + (h*J[[1,2]])/J[[1,1]], rhat-h, mhat - (h*J[[1,2]])/J[[1,1]], rhat+h}] (* -log of likelihood function using the Lea-Coulson normal approximation *) leaLike[m_, data_] := Plus@@((11.6/(#1/m - Log[m] + 4.5) - 2.02)&/@data) (* ==================== zero counterparts ======================= *) (* zero counterpart of Luria-Delbruck distribution *) pmfLD0[m_/;Positive[m], phi_/;(phi>0)&&(phi<=1), k_Integer/;NonNegative[k]] := pmfLuria0[N[m],N[phi],k] pmfLD0[m_, k_Integer]:=pmfLuria0[m,1,k]/;(m>0)&&(k>=0) (* zero counterpart of using Newton-Raphson method *) newtonLD0[data_List, opts___] := Module[ {m0, m1, score, info, showIter, phi, tol, maxIter, initEst, i}, showIter = ShowIterations /. {opts} /. Options[newtonLD0]; initEst = InitialM /. {opts} /. Options[newtonLD0]; maxIter = MaxIterations /. {opts} /. Options[newtonLD0]; phi = Phi /. {opts} /. Options[newtonLD0]; tol = Tolerance /. {opts} /. Options[newtonLD0]; If[initEst>0, m0=initEst, If[Median[data]==0, m0=MethodOfP0[data], m0=jonesMedian[data]]]; If[showIter, Print[m0]]; For[i = 1, i <= maxIter, i++, {score, info} = LDScoreAndInfo0[m0, phi, data]; m1 = m0 + score/info; If[Abs[(m1 - m0)/m0] < tol, Return[m1], If[showIter, Print[m1]]; m0 = m1]]; Print["no convergence!"]; ] (* zero counterpart of asymptotic CI for m with LD distribution, 5/4/04, 4/24/06 *) LRIntervalLD0[data_List, opts___] := Module[ {mhat, la, qa, h, m0, m1, like, score, info, mL, mU, alpha, tol, initM, initLow, initUp, maxIter, phi, showIter}, alpha = Alpha /. {opts} /. Options[LRIntervalLD0]; tol = Tolerance /. {opts} /. Options[LRIntervalLD0]; maxIter = MaxIterations /. {opts} /. Options[LRIntervalLD0]; phi = Phi /. {opts} /. Options[LRIntervalLD0]; initM= InitialM /. {opts} /. Options[LRIntervalLD0]; initLow= InitialLowerM /. {opts} /. Options[LRIntervalLD0]; initUp= InitialUpperM /. {opts} /. Options[LRIntervalLD0]; showIter= ShowIterations /. {opts} /. Options[LRIntervalLD0]; If[phi<0||phi>1, Return["\!\(\[Phi]\) out of range!"]]; mhat = newtonLD0[data, Phi->phi, InitialM->initM, Tolerance->tol]; {like, score} = LDLoglikeScore0[mhat, phi, data]; {score, info} =LDScoreAndInfo0[mhat,phi,data]; qa=Quantile[ChiSquareDistribution[1], 1-alpha]; h=Sqrt[qa/info]; la = like - 0.5*qa; If[showIter, Print["Iterating for lower limit ..."]]; If[initLow<=0, m0=mhat-0.5*h, m0=initLow]; For[i = 1, i <= maxIter, i++, {like, score} = LDLoglikeScore0[m0, phi, data]; m1 = m0 - (like - la)/score; If[showIter, Print[m0]]; If[Abs[(m1 - m0)/m0] <= 10^(-5), mL=m1; Break[], m0 = m1]]; If[showIter, Print["Iterating for upper limit ..."]]; If[initUp<=0, m0=mhat+0.5*h, m0=initUp]; For[i = 1, i <= maxIter, i++, {like, score} = LDLoglikeScore0[m0, phi, data]; m1 = m0 - (like - la)/score; If[showIter, Print[m0]]; If[Abs[(m1 - m0)/m0] <= tol, mU=m1; Break[], m0 = m1]]; {mL, mU}] (* ============== internal functions of zero counterparts ============== *) (* -------- probability mass function of LD dist ---------- *) pmfLuria0[m_, phi_, k_Integer] := Module[{cumuphi, prob, n, j}, prob = Table[0, {k + 1}]; prob[[1]] = Exp[-N[m]]; For[n = 1, n <= k, n++, cumuphi = 1.; For[j = 1, j <= n, j++, prob[[n + 1]] += cumuphi*(1. - j*(phi/(j + 1)))*prob[[n - j + 1]]; cumuphi *= phi; ]; prob[[n + 1]] = prob[[n + 1]]*(m/n); ]; prob] (* ------------- sequence convolution ------------- *) pdfConv0[x_List, y_List] := Module[{nx, z, k, i}, nx = Length[x]; z = Table[0, {nx}]; For[k = 1, k <= nx, k++, For[i = 1, i <= k, i++, z[[k]] += x[[i]]*y[[k - i + 1]]; ]]; z] /; Length[x] == Length[y] (* --------- Fisher's information ------------ *) LDScoreAndInfo0[m_?NumericQ, phi_, data_List] := Module[{n, h, p, p1, p2, info, score}, n = Max[data]; p = pmfLuria0[N[m], N[phi], n]; h = Table[phi^(k - 1)*(1/k - phi/(k + 1)), {k, 1, n}]; h = Prepend[h, -1]; p1 = pdfConv0[p, N[h]]; p2 = pdfConv0[p1, N[h]]; score = ((p1/p)[[#1 + 1]] & ) /@ data; info = (((p1/p)^2 - p2/p)[[#1 + 1]] & ) /@ data; {Plus @@ score, Plus @@ info}] (* loglikelihood and score, 5-3-04, 4-24-06 *) LDLoglikeScore0[m_?NumericQ,phi_,data_List] := Module[ {n, h, p, p1, loglike, score}, n= Max[data]; p = pmfLuria0[N[m], N[phi], n]; loglike = Log[p[[#1 + 1]]] & /@ data; h = Table[phi^(k-1)*(1/k - phi/(k + 1)), {k, 1, n}]; h = Prepend[h, -1]; p1 = pdfConv0[p, N[h]]; score = ((p1/p)[[#1 + 1]])& /@ data; {Plus @@ loglike, Plus @@ score} ] (* ================== end of zero counterparts ================== *) (* =================== functions for Haldane's formulation ================ *) (* numbers of mutations under Haldane's formulation *) haldaneMutations[mu_, g_, N0_:1] := mu*N0*(((2 - mu)^g - 1)/(1 - mu)) (* simulation of synchronized mutation model *) simuHaldane[gen_Integer, mu_, opts___] := Module[ {wild = 1, mut = 0, newWild, newMut, growth}, growth = ShowGrowth /. {opts} /. Options[simuHaldane]; For[i = 1, i <= gen, i++, newMut = newWild = 0; For[j = 1, j <= wild, j++, If[Random[] < mu, newMut++, newWild++]]; mut = 2*mut + newMut; wild = wild + newWild; If[growth, Print["gen=", i, " wild=", wild, " mutant=", mut] ]; ]; mut] (* Angerer's simulation scheme *) simuAngerer[mu_, div_, opts___] := Module[{wild, mut = 0, who, init, growth}, wild = InitialSize /. {opts} /. Options[simuAngerer]; growth = ShowGrowth /. {opts} /. Options[simuAngerer]; For[i = 1, i <= div, i++, who = mut/(mut + wild); If[Random[] < who, mut++, If[Random[] < mu, mut++, wild++]]; If[growth, Print["div=", i, " wild=", wild, " mutant=", mut]]]; mut] (* ********* Haldane's formulation, mathlink part added January 20, 2008 *******) (* -------- slow, Haldane ----------- *) pmfHaldane0[gen_, mu_, n_, N0_:1] := Module[ {p, prob, A, nn, g, k, j}, nn = Reverse[NestList[Floor[#1/2] & , n, gen - 1]]; p = prob = Table[0, {n + 1}]; prob[[Range[nn[[1]]+1]]] = Table[PDF[BinomialDistribution[N0,mu],j],{j, 0, nn[[1]]}]; For[g = 2, g <= gen, g++, p = prob; prob = Table[0, {n + 1}]; For[k = 0, k <= nn[[g]], k++, A = getAkj0[g - 1, k, mu, N0]; For[j = Max[0, k - 2^(g - 1)*N0], j <= Floor[k/2], j++, prob[[k + 1]] = prob[[k + 1]] + A[[j + 1]]*p[[j + 1]]]]]; prob] getAkj0[g_, k_, mu_, N0_] := Module[{Ng = 2^g*N0, U = Floor[k/2]}, Akj = Table[0, {U + 1}]; Akj[[U + 1]] = AkHalfk[g, mu, k, N0]; For[j = U - 1, j >= 0, j--, Akj[[j + 1]] = (Ng - j)*(Ng - k + j + 1)* (mu^2/((k - 2*j)*(k - 2*j - 1)*(1 - mu)))*Akj[[j + 2]]; ]; Akj] AkHalfk[g_, mu_, k_, N0_] := Module[{Ng}, Ng = N0*2^g; If[EvenQ[k], (1 - mu)^(Ng - k/2), (Ng - (k - 1)/2)*mu*(1 - mu)^(Ng - (k + 1)/2)]] (* =========== MLE for Haldane ========== *) newtonHaldane0[data_, g_, opts___] := Module[{mu0, mu1, U, J, i, tol, showIter, maxIter, initEst, N0}, showIter = ShowIterations /. {opts} /. Options[newtonHaldane0]; initEst = InitialMu /. {opts} /. Options[newtonHaldane0]; tol = Tolerance /. {opts} /. Options[newtonHaldane0]; maxIter = MaxIterations /. {opts} /. Options[newtonHaldane0]; N0 = InitialCells /. {opts} /. Options[newtonHaldane0]; mu0 = If[initEst > 0, initEst, 0.5*rossman[data, g, N0]]; If[showIter, Print[mu0]]; For[i = 1, i <= maxIter, i++, {U, J} = getUAndJ0[data, g, mu0, N0]; mu1 = mu0 + U/J; If[Abs[(mu1 - mu0)/mu0] < tol, Return[mu1], If[showIter, Print[mu1]]; mu0 = mu1; ]]; Print["no convergence!"]; ] (* ------------ auxiliary functions for newtonHaldane ---------- *) (* ----------- the a sequence and its first derivative ---------- *) a[g_, k_, j_, mu_, N0_] := (k - 2*j - (N0*2^g - j)*mu)/(mu*(1 - mu)) ap[g_, k_, j_, mu_, N0_] := (2*j - k + 2*(k - 2*j)*mu - (N0*2^g - j)*mu^2)/ (mu^2*(1 - mu)^2) (* ---- computing first and second derivatives of the Haldane dist --- *) derivHaldane0[gen_, mu_, n_, N0_] := Module[ {prob0, prob10, prob20, prob, prob1, prob2, A, A0, Ng = 2^gen*N0, nn, j, k}, nn = Reverse[NestList[Floor[#1/2] & , n, gen - 1]]; prob0 = prob = prob10 = prob1 = prob20 = prob2 = Table[0, {n + 1}]; prob[[Range[nn[[1]] + 1]]] = Table[PDF[BinomialDistribution[N0,mu],j], {j, 0, nn[[1]]}]; prob1[[Range[nn[[1]] + 1]]] = Table[j - N0*mu, {j, 0, nn[[1]]}]*prob[[Range[nn[[1]] + 1]]]; prob1 = prob1/(mu*(1 - mu)); prob2[[Range[nn[[1]] + 1]]] = Table[j*(j - 1) - 2*j*(N0 - 1)*mu + N0*(N0 - 1)*mu^2, {j, 0, nn[[1]]}]* prob[[Range[nn[[1]] + 1]]]; prob2 = prob2/(mu*(1 - mu))^2; For[g = 2, g <= gen, g++, prob0 = prob; prob10 = prob1; prob20 = prob2; prob = prob1 = prob2 = Table[0, {n + 1}]; For[k = 0, k <= nn[[g]], k++, A = getAkj0[g - 1, k, mu, N0]; For[j = Max[0, k - 2^(g - 1)*N0], j <= Floor[k/2], j++, prob[[k + 1]] += A[[j + 1]]*prob0[[j + 1]]; prob1[[k + 1]] += A[[j + 1]]*a[g - 1, k, j, mu, N0]*prob0[[ j + 1]] + A[[j + 1]]*prob10[[j + 1]]; prob2[[k + 1]] += A[[j + 1]]*(a[g - 1, k, j, mu, N0]^2 + ap[g - 1, k, j, mu, N0])* prob0[[j + 1]] + 2*A[[j + 1]]*a[g - 1, k, j, mu, N0]* prob10[[ j + 1]] + A[[j + 1]]*prob20[[j + 1]]]]]; {prob, prob1, prob2}] (* ------------ computing score and Fisher's information ----------------- *) getUAndJ0[data_List, g_, mu_, N0_] := Module[{n, p, p1, p2, U, J}, n = Max[data]; {p, p1, p2} = derivHaldane0[g, mu, n, N0]; p = p[[data + 1]]; p1 = p1[[data + 1]]; p2 = p2[[data + 1]]; U = Plus @@ (p1/p); J = Plus @@ ((p1/p)^2 - p2/p); {U, J}] (* -------- the method of the mean, by Rossman, for initial guess ---------- *) rossman[data_List, g_Integer, N0_Integer:1] := Module[{m, Nt}, Nt = N0*2^g; m = Mean[data]; 2*(1. - (1 - m/Nt)^(1/g))] (* ============= confidence interval for Haldane's dist ========== *) CIHaldane0[data_, g_, opts___] := Module[ {N0, alpha, maxIter, tol, showIter, initmu, initLow, initUp, muHat, h, chiAlphal, lik, U, J, la, ci0, ci1, i, muLow, muUp}, N0 = InitialCells /. {opts} /. Options[CIHaldane0]; alpha = Alpha /. {opts} /. Options[CIHaldane0]; maxIter = MaxIterations /. {opts} /. Options[CIHaldane0]; tol = Tolerance /. {opts} /. Options[CIHaldane0]; showIter = ShowIterations /. {opts} /. Options[CIHaldane0]; initmu = InitialMu /. {opts} /. Options[CIHaldane0]; initLow = InitialLowerLimit /. {opts} /. Options[CIHaldane0]; initUp = InitialUpperLimit /. {opts} /. Options[CIHaldane0]; chiAlpha = Quantile[ChiSquareDistribution[1], 1 - alpha]; If[showIter, Print["Iterating for MLE of \!\(\[Mu]\)..."]]; If[initmu < 0, muHat = newtonHaldane0[data, g, InitialCells -> N0, ShowIterations->showIter], muHat = newtonHaldane0[data, g, InitialCells -> N0, InitialMu -> initmu, ShowIterations->showIter]]; {U, J} = getUAndJ0[data, g, muHat, N0]; h = 0.5*Sqrt[chiAlpha/J]; {lik, U} = getLikAndScore0[data, g, muHat, N0]; la = lik - chiAlpha/2; If[showIter, Print["Iterating for lower limit of CI ..."]]; If[initLow < 0, ci0 = muHat - h, ci0 = initLow]; For[i = 1, i <= maxIter, i++, {lik, U} = getLikAndScore0[data, g, ci0, N0]; ci1 = ci0 - (lik - la)/U; If[showIter, Print[ci0]]; If[Abs[(ci1 - ci0)/ci0] < tol, muLow = ci1; Break[], ci0 = ci1]; ]; If[showIter, Print["Iterating for upper limit of CI ..."]]; If[initUp < 0, ci0 = muHat + h, ci0 = initUp]; For[i = 1, i <= maxIter, i++, {lik, U} = getLikAndScore0[data, g, ci0, N0]; ci1 = ci0 - (lik - la)/U; If[showIter, Print[ci0]]; If[Abs[(ci1 - ci0)/ci0] < tol, muUp = ci1; Break[], ci0 = ci1]; ]; {muLow, muUp}] getUAndJ0[data_List, g_, mu_, N0_] := Module[{n, p, p1, p2, U, J}, n = Max[data]; {p, p1, p2} = derivHaldane0[g, mu, n, N0]; p = p[[data + 1]]; p1 = p1[[data + 1]]; p2 = p2[[data + 1]]; U = Plus @@ (p1/p); J = Plus @@ ((p1/p)^2 - p2/p); {U, J}] getLikAndScore0[data_List, g_, mu_, N0_] := Module[{n, p, p1, U, lik}, n = Max[data]; {p, p1} = pAndP1Haldane0[g, mu, n, N0]; p = p[[data + 1]]; p1 = p1[[data + 1]]; lik = Plus @@ Log[p]; U = Plus @@ (p1/p); {lik, U}] pAndP1Haldane0[gen_, mu_, n_, N0_:1] := Module[{prob0, prob10, prob, prob1, A, A0, Ng = 2^gen*N0, nn, j, k}, nn = Reverse[NestList[Floor[#1/2] & , n, gen - 1]]; prob0 = prob = prob10 = prob1 = Table[0, {n + 1}]; prob[[Range[nn[[1]] + 1]]] = (* typo corrected 4:47 pm, Jan 22, 2008 *) Table[PDF[BinomialDistribution[N0, mu],j], {j, 0, nn[[1]]}]; prob1[[Range[nn[[1]] + 1]]] = Table[j - N0*mu, {j, 0, nn[[1]]}]*prob[[Range[nn[[1]] + 1]]]; prob1 = prob1/(mu*(1 - mu)); For[g = 2, g <= gen, g++, prob0 = prob; prob10 = prob1; prob = prob1 = Table[0, {n + 1}]; For[k = 0, k <= nn[[g]], k++, A = getAkj0[g - 1, k, mu, N0]; For[j = Max[0, k - 2^(g - 1)*N0], j <= Floor[k/2], j++, prob[[k + 1]] += A[[j + 1]]*prob0[[j + 1]]; prob1[[k + 1]] += A[[j + 1]]*a[g - 1, k, j, mu, N0]*prob0[[j + 1]] + A[[j + 1]]*prob10[[j + 1]]; ]]]; {prob, prob1}] (* ------------- end of CI for Haldane ------------------- *) (* =========== Drake's formula under Haldane's model, May 21, 2006 ===== *) drakeHaldane[data_,Ng_, opts___] := Module[ {med, m0, m1, showIter, guess, tol, maxIter, useMed, i}, showIter = ShowIterations /. {opts} /. Options[drakeHaldane]; guess = initialM /. {opts} /. Options[drakeHaldane]; tol = Tolerance /. {opts} /. Options[drakeHaldane]; maxIter = MaxIterations /. {opts} /. Options[drakeHaldane]; useMed = UseMedian /. {opts} /. Options[drakeHaldane]; med =If[useMed, Median[data], Mean[data]]*Log[2]*2; If[guess<=0, m0=1.0, m0 = Ng*N[guess]]; If[showIter, Print[m0/Ng]]; For[i = 1, i <= maxIter, i++, m1 = m0 - (m0*Log[m0]-med)/(1+Log[m0]); If[showIter, Print[m1/Ng]]; If[Abs[(m0 - m1)/m0] < tol, Return[m1/Ng], m0 = m1]]; Print["no convergence"];] (* -------- end of Drake's formula ----------- *) (* ----------- Haldane, mathlink, fast ------------- *) pmfHaldane[gen_Integer, mu_, n_Integer, N0_Integer:1]:= $luria`$pmfHaldane[gen,mu,n,N0]/;(gen>=1)&&(mu>0)&&(n>=0)&&(n<=N0*2^gen) (* =========== MLE for Haldane ========== *) newtonHaldane[data_, g_, opts___] := Module[{mu0, mu1, U, J, i, tol, showIter, maxIter, initEst, N0}, showIter = ShowIterations /. {opts} /. Options[newtonHaldane]; initEst = InitialMu /. {opts} /. Options[newtonHaldane]; tol = Tolerance /. {opts} /. Options[newtonHaldane]; maxIter = MaxIterations /. {opts} /. Options[newtonHaldane]; N0 = InitialCells /. {opts} /. Options[newtonHaldane]; mu0 = If[initEst > 0, initEst, 0.5*rossman[data, g, N0]]; If[showIter, Print[mu0]]; For[i = 1, i <= maxIter, i++, {U, J} = getUAndJ[data, g, mu0, N0]; mu1 = mu0 + U/J; If[Abs[(mu1 - mu0)/mu0] < tol, Return[mu1], If[showIter, Print[mu1]]; mu0 = mu1; ]]; Print["no convergence!"]; ] (* ------------ auxiliary functions for newtonHaldane ---------- *) (* ------------ computing score and Fisher's information ----------------- *) getUAndJ[data_List, g_, mu_, N0_] := Module[{n, p, p1, p2, U, J}, n = Max[data]; {p, p1, p2} = $luria`$derivHaldane[g, mu, n, N0]; p = p[[data + 1]]; p1 = p1[[data + 1]]; p2 = p2[[data + 1]]; U = Plus @@ (p1/p); J = Plus @@ ((p1/p)^2 - p2/p); {U, J}] (* ============= confidence interval for Haldane's dist ========== *) CIHaldane[data_, g_, opts___] := Module[ {N0, alpha, maxIter, tol, showIter, initmu, initLow, initUp, muHat, h, chiAlphal, lik, U, J, la, ci0, ci1, i, muLow, muUp}, N0 = InitialCells /. {opts} /. Options[CIHaldane]; alpha = Alpha /. {opts} /. Options[CIHaldane]; maxIter = MaxIterations /. {opts} /. Options[CIHaldane]; tol = Tolerance /. {opts} /. Options[CIHaldane]; showIter = ShowIterations /. {opts} /. Options[CIHaldane]; initmu = InitialMu /. {opts} /. Options[CIHaldane]; initLow = InitialLowerLimit /. {opts} /. Options[CIHaldane]; initUp = InitialUpperLimit /. {opts} /. Options[CIHaldane]; chiAlpha = Quantile[ChiSquareDistribution[1], 1 - alpha]; If[showIter, Print["Iterating for MLE of \!\(\[Mu]\)..."]]; If[initmu < 0, muHat = newtonHaldane[data, g, InitialCells -> N0, ShowIterations->showIter], muHat = newtonHaldane[data, g, InitialCells -> N0, InitialMu -> initmu, ShowIterations->showIter]]; (* If[showIter, Print["mle for \!\(\[Mu]\) is ... ", muHat]]; *) {U, J} = getUAndJ[data, g, muHat, N0]; h = 0.5*Sqrt[chiAlpha/J]; {lik, U} = getLikAndScore[data, g, muHat, N0]; la = lik - chiAlpha/2; If[showIter, Print["Iterating for lower limit of CI ..."]]; If[initLow < 0, ci0 = muHat - h, ci0 = initLow]; For[i = 1, i <= maxIter, i++, {lik, U} = getLikAndScore[data, g, ci0, N0]; ci1 = ci0 - (lik - la)/U; If[showIter, Print[ci0]]; If[Abs[(ci1 - ci0)/ci0] < tol, muLow = ci1; Break[], ci0 = ci1]; ]; If[showIter, Print["Iterating for upper limit of CI ..."]]; If[initUp < 0, ci0 = muHat + h, ci0 = initUp]; For[i = 1, i <= maxIter, i++, {lik, U} = getLikAndScore[data, g, ci0, N0]; ci1 = ci0 - (lik - la)/U; If[showIter, Print[ci0]]; If[Abs[(ci1 - ci0)/ci0] < tol, muUp = ci1; Break[], ci0 = ci1]; ]; {muLow, muUp}] getLikAndScore[data_List, g_, mu_, N0_] := Module[{n, p, p1, U, lik}, n = Max[data]; {p, p1} = $luria`$pAndP1Haldane[g, mu, n, N0]; p = p[[data + 1]]; p1 = p1[[data + 1]]; lik = Plus @@ Log[p]; U = Plus @@ (p1/p); {lik, U}] (* ============== end of functions for Haldane's formulation ============ *) (* **************** Bartlett's formulation ******************* *) (* --------------- simulation ---------------- *) simuBartlett[b_, mu_, N0_, T_, opts___] := Module[{clock, wild, mutant, arrival, count, events, showEve, maxEve}, maxEve = MaxEvents /. {opts} /. Options[simuBartlett]; showEve = ShowGrowth /. {opts} /. Options[simuBartlett]; wild = N0; mutant = 0; clock = 0; events = 0; While[True, If[showEve, Print[StringJoin["wild= ", ToString[wild], " mutant=", ToString[mutant]]]]; arrival = -Log[Random[]]/((wild + mutant)*b); If[clock + arrival > T, Return[{wild, mutant}], clock += arrival]; If[Random[] > wild/(wild + mutant), mutant += 1, If[Random[] > mu/b, wild += 1, mutant += 1]]; events += 1; If[events > maxEve, Print["limit reached, simulation abandoned."]; Return[]]]; ] (* -------------- a new distribution, Nov 29, 2007 --------------- *) pmfBli[A_, n_] := Module[{q, prob}, prob = q = Table[0, {n + 1}]; q[[1]] = 1 + A; For[i = 1, i <= n, i++, q[[i + 1]] = (-A)*(1/i - 1/(i + 1))]; prob[[1]] = 1/(1 + A); For[i = 1, i <= n, i++, For[j = 0, j <= i - 1, j++, prob[[i + 1]] += prob[[j + 1]]*q[[i - j + 1]]]; prob[[i + 1]] /= -q[[1]]]; prob] (* ---------- mean function ----------- *) bartlettMean[x_, N0_, phi_] := Module[{xbar}, xbar=If[Median[x]>0, N[Median[x]], N[Mean[x]] ]; Log[1 - (1-phi)*xbar/N0]/Log[1-phi] ] (* ------------ pmf ---------------- *) pmfBartlett[a_,phi_,n_]:=$luria`$pmfBartlett[a,phi,n] pmfBartlett[a_,phi_,n_,k_Integer/;k>0]:=Module[{p}, p=$luria`$pmfBartlett[a,phi,n]; Nest[$luria`$pdfconv[p,#]&,p,k-1] ] (* ----- Newton-Raphson for finding mle of mutation rates, May 6, 2007 -------- *) newtonBartlett[data_List, N0_, phi_, opts___] := Module[{a0, a1, score, info, showIter, tol, maxIter, initGuess,i}, initGuess = InitialMutationRate /. {opts} /. Options[newtonBartlett]; showIter = ShowIterations /. {opts} /. Options[newtonBartlett]; maxIter = MaxIterations /. {opts} /. Options[newtonBartlett]; tol = Tolerance /. {opts} /. Options[newtonBartlett]; If[initGuess > 0, a0 = N[initGuess], a0 = bartlettMean[data, N0, phi]]; If[showIter, Print[a0]]; For[i = 1, i <= maxIter, i++, {score, info} = scoreInfoBart[a0, phi, N0, data]; a1 = a0 + score/info; If[Abs[(a1 - a0)/a0] < tol, Return[a1], If[showIter, Print[a1]]; a0 = a1; ]]; Print["no convergence!"]; ] scoreInfoBart[alpha_, phi_, N0_, data_List] := Module[{n, h, p, p1, p2, score, info}, n = Max[data]; p = pmfBartlett[alpha, phi, n, N0]; {p1, p2} = getDerivativeBart[alpha, phi, N0, n]; score = ((p1/p)[[#1 + 1]] & ) /@ data; info = (((p1/p)^2 - p2/p)[[#1 + 1]] & ) /@ data; {Plus @@ score, Plus @@ info}] getDerivativeBart[alpha_, phi_, N0_, n_] := Module[{A, p0, p1, p2, pN, c, k}, A = getAk[alpha, phi, N0, n]; p0 = pmfBartlett[alpha, phi, n]; c = Table[phi^k/k, {k, 1, n}]; c = Prepend[c, 0]; c = ((N0 + 1)/(1 - phi))*$luria`$pdfconv[N[A], p0] - c; pN = pmfBartlett[alpha, phi, n, N0 + 1]; p1 = $luria`$pdfconv[N[A], pN]*(N0/(1 - phi)); p2 = $luria`$pdfconv[N[p1], N[c]]; {p1, p2}] getAk[alpha_, phi_, N0_, K_] := Module[{a, b, p, pp, i, allphi}, allphi = Table[phi^i, {i, 0, K}]; a = Table[1/i - phi/(i + 1), {i, 1, K}]; a = Prepend[a, -phi]; b = Table[0, {K + 1}]; b[[1]] = 1; For[i = 1, i <= K, i++, b[[i + 1]] = ((i - 1 - alpha)/i)*b[[i]]; ]; allphi*$luria`$pdfconv[N[a], N[b]] ] (* ------------ Newton-Raphson for CI under Bartlett formulation ------------ *) CIBartlett[data_List, N0_, phi_, opts___] := Module[{mhat, la, qa, h, m0, m1, like, score, info, mLow, mUp, alpha, tol, initMu, initLow, initUp, maxIter, showIter, i}, alpha = Alpha /. {opts} /. Options[CIBartlett]; tol = Tolerance /. {opts} /. Options[CIBartlett]; showIter = ShowIterations /. {opts} /. Options[CIBartlett]; maxIter = MaxIterations /. {opts} /. Options[CIBartlett]; initMu = InitialMutationRate /. {opts} /. Options[CIBartlett]; initLow = InitialLowerLimit /. {opts} /. Options[CIBartlett]; initUp = InitialUpperLimit /. {opts} /. Options[CIBartlett]; If[showIter, Print["Iterating for MLE ..."]]; mhat = newtonBartlett[data, N0, phi, InitialMutationRate -> initMu, ShowIterations -> showIter]; {like, score} = loglikScoreBart[mhat, phi, N0, data]; {score, info} = scoreInfoBart[mhat, phi, N0, data]; qa = Quantile[ChiSquareDistribution[1], 1 - alpha]; h = Sqrt[qa/info]; la = like - 0.5*qa; If[showIter, Print["Iterating for lower limit ..."]]; If[initLow <= 0, m0 = mhat - 0.5*h, m0 = initLow]; For[i = 1, i <= maxIter, i++, {like, score} = loglikScoreBart[m0, phi, N0, data]; m1 = m0 - (like - la)/score; If[showIter, Print[m0]]; If[Abs[(m1 - m0)/m0] < tol, mLow = m1; Break[], m0 = m1]]; If[showIter, Print["Iterating for upper limit ..."]]; If[initUp <= 0, m0 = mhat + 0.5*h, m0 = initUp]; For[i = 1, i <= maxIter, i++, {like, score} = loglikScoreBart[m0, phi, N0, data]; m1 = m0 - (like - la)/score; If[showIter, Print[m0]]; If[Abs[(m1 - m0)/m0] < tol, mUp = m1; Break[], m0 = m1]]; {mLow, mUp}] loglikScoreBart[alpha_, phi_, N0_, data_List] := Module[{n, h, p, p1, p2, score, loglik}, n = Max[data]; p = pmfBartlett[alpha, phi, n, N0]; loglik = (Log[p[[#1 + 1]]] & ) /@ data; {p1, p2} = getDerivativeBart[alpha, phi, N0, n]; score = ((p1/p)[[#1 + 1]] & ) /@ data; {Plus @@ loglik, Plus @@ score}] (* ------ end of Newton-Raphson ---------- *) (* ----------- golden section search for MLE, April 20, 2007 ------- *) goldenBartlett[x_, N0_, phi_, opts___] := Module[{lik, mom, maxIter, initTrip, tol, showIter}, lik = Function[{z}, -likely[x, z, phi, N0]]; initTrip = InitialTriplet /. {opts} /. Options[goldenBartlett]; If[initTrip == {1, 1, 1}, mom = bartlettMean[x, N0, phi]; initTrip = {mom/10, mom, 10*mom}]; initTrip=N[initTrip]; If[!( (lik[initTrip[[1]]]>lik[initTrip[[2]]])&& (lik[initTrip[[3]]]>lik[initTrip[[2]]]) ), Print["Invalid initial triplet ..."]; Return[] ]; maxIter = MaxIterations /. {opts} /. Options[goldenBartlett]; tol = Tolerance /. {opts} /. Options[goldenBartlett]; showIter = ShowIterations /. {opts} /. Options[goldenBartlett]; golden[lik, initTrip, maxIter, tol, showIter]] likely[x_, a_, phi_, N0_] := Module[{p, n = Max[x]}, p = pmfBartlett[a, phi, n, N0]; Plus @@ Log[p[[x + 1]]]] golden[f_, {a0_, b0_, c0_}, iter_, tol_, showIter_] := Module[{fa, fb, fc, f0, x0, a = a0, b = b0, c = c0, i}, fa = f[a]; fb = f[b]; fc = f[c]; For[i = 1, i <= iter, i++, x0 = If[c - b > b - a, b + 0.381966*(c - b), b - 0.381966*(b - a)]; f0 = f[x0]; If[f0 > fb && x0 > b, {a, b, c} = {a, b, x0}; fc = f0, If[f0 < fb && x0 > b, {a, b, c} = {b, x0, c}; {fa, fb, fc} = {fb, f0, fc}, If[f0 > fb && x0 < b, {a, b, c} = {x0, b, c}; {fa, fb, fc} = {f0, fb, fc}, {a, b, c} = {a, x0, b}; {fa, fb, fc} = {fa, f0, fb}]]]; If[showIter, Print[StringJoin["i=", ToString[i], ", a=", ToString[TraditionalForm[a]], ", b=", ToString[TraditionalForm[ b]], ", c=", ToString[TraditionalForm[c]], ", -log(L)=", ToString[fb]]]]; If[Abs[c - a]/Abs[b] < tol, Return[b]]; ]; Print["Golden section search fails to converge ..."]; ] (* --------- zero (0) version, or the slow version, April 26, 2007 ------- *) (* --------- zero (0) version, modified January 21, 2008 ------- *) listConv[x_List, y_List] := Module[{n = Length[x]}, (x[[Range[#1]]] . Reverse[y[[Range[#1]]]] & ) /@ Range[n]] pmfBartlett0[a_,phi_,n_,k_Integer/;k>0]:=Module[{p}, p=pmfBartlett0[a,phi,n]; Nest[listConv[p,#]&,p,k-1] ] pmfBartlett0[a_, phi_, n_] := Module[{b, t, p, k, j}, b = Table[0, {n}]; b[[1]] = (-a)*phi; For[k = 1, k < n, k++, b[[k + 1]] = phi*((k - a)/(k + 1))*b[[k]]]; t = Table[0, {n}]; For[k = 1, k <= n, k++, t[[k]] = b[[k]]*(phi*((a - k)/(k + 1)) + 1)]; t = Prepend[t, 1 - (1 - a)*phi]; p = Table[0, {n + 1}]; p[[1]] = (1 - phi)/(1 - (1 - a)*phi); For[k = 1, k <= n, k++, p[[k + 1]] = 0; For[j = 0, j <= k - 1, j++, p[[k + 1]] = p[[k + 1]] + p[[j + 1]]*t[[k - j + 1]]]; p[[k + 1]] = -p[[k + 1]]/t[[1]]; ]; p] (* ------------- *) newtonBartlett0[data_List, N0_, phi_, opts___] := Module[{a0, a1, score, info, showIter, tol, maxIter, initGuess,i}, initGuess = InitialMutationRate /. {opts} /. Options[newtonBartlett0]; showIter = ShowIterations /. {opts} /. Options[newtonBartlett0]; maxIter = MaxIterations /. {opts} /. Options[newtonBartlett0]; tol = Tolerance /. {opts} /. Options[newtonBartlett0]; If[initGuess > 0, a0 = N[initGuess], a0 = bartlettMean[data, N0, phi]]; If[showIter, Print[a0]]; For[i = 1, i <= maxIter, i++, {score, info} = scoreInfoBart0[a0, phi, N0, data]; a1 = a0 + score/info; If[Abs[(a1 - a0)/a0] < tol, Return[a1], If[showIter, Print[a1]]; a0 = a1; ]]; Print["no convergence!"]; ] scoreInfoBart0[alpha_, phi_, N0_, data_List] := Module[{n, h, p, p1, p2, score, info}, n = Max[data]; p = pmfBartlett0[alpha, phi, n, N0]; {p1, p2} = getDerivativeBart0[alpha, phi, N0, n]; score = ((p1/p)[[#1 + 1]] & ) /@ data; info = (((p1/p)^2 - p2/p)[[#1 + 1]] & ) /@ data; {Plus @@ score, Plus @@ info}] getDerivativeBart0[alpha_, phi_, N0_, n_] := Module[{A, p0, p1, p2, pN, c, k}, A = getAk0[alpha, phi, N0, n]; p0 = pmfBartlett0[alpha, phi, n]; c = Table[phi^k/k, {k, 1, n}]; c = Prepend[c, 0]; c = ((N0 + 1)/(1 - phi))* listConv[N[A], p0] - c; pN = pmfBartlett0[alpha, phi, n, N0 + 1]; p1 = listConv[N[A], pN]*(N0/(1 - phi)); p2 = listConv[N[p1], N[c]]; {p1, p2}] getAk0[alpha_, phi_, N0_, K_] := Module[{a, b, p, pp, i, allphi}, allphi = Table[phi^i, {i, 0, K}]; a = Table[1/i - phi/(i + 1), {i, 1, K}]; a = Prepend[a, -phi]; b = Table[0, {K + 1}]; b[[1]] = 1; For[i = 1, i <= K, i++, b[[i + 1]] = ((i - 1 - alpha)/i)*b[[i]]; ]; allphi* listConv[N[a], N[b]] ] (* ------------ Newton-Raphson for CI under Bartlett formulation ------------ *) CIBartlett0[data_List, N0_, phi_, opts___] := Module[{mhat, la, qa, h, m0, m1, like, score, info, mLow, mUp, alpha, tol, initMu, initLow, initUp, maxIter, showIter, i}, alpha = Alpha /. {opts} /. Options[CIBartlett0]; tol = Tolerance /. {opts} /. Options[CIBartlett0]; showIter = ShowIterations /. {opts} /. Options[CIBartlett0]; maxIter = MaxIterations /. {opts} /. Options[CIBartlett0]; initMu = InitialMutationRate /. {opts} /. Options[CIBartlett0]; initLow = InitialLowerLimit /. {opts} /. Options[CIBartlett0]; initUp = InitialUpperLimit /. {opts} /. Options[CIBartlett0]; If[showIter, Print["Iterating for MLE ..."]]; mhat = newtonBartlett0[data, N0, phi, InitialMutationRate -> initMu, ShowIterations -> showIter]; {like, score} = loglikScoreBart0[mhat, phi, N0, data]; {score, info} = scoreInfoBart0[mhat, phi, N0, data]; qa = Quantile[ChiSquareDistribution[1], 1 - alpha]; h = Sqrt[qa/info]; la = like - 0.5*qa; If[showIter, Print["Iterating for lower limit ..."]]; If[initLow <= 0, m0 = mhat - 0.5*h, m0 = initLow]; For[i = 1, i <= maxIter, i++, {like, score} = loglikScoreBart0[m0, phi, N0, data]; m1 = m0 - (like - la)/score; If[showIter, Print[m0]]; If[Abs[(m1 - m0)/m0] < tol, mLow = m1; Break[], m0 = m1]]; If[showIter, Print["Iterating for upper limit ..."]]; If[initUp <= 0, m0 = mhat + 0.5*h, m0 = initUp]; For[i = 1, i <= maxIter, i++, {like, score} = loglikScoreBart0[m0, phi, N0, data]; m1 = m0 - (like - la)/score; If[showIter, Print[m0]]; If[Abs[(m1 - m0)/m0] < tol, mUp = m1; Break[], m0 = m1]]; {mLow, mUp}] loglikScoreBart0[alpha_, phi_, N0_, data_List] := Module[{n, h, p, p1, p2, score, loglik}, n = Max[data]; p = pmfBartlett0[alpha, phi, n, N0]; loglik = (Log[p[[#1 + 1]]] & ) /@ data; {p1, p2} = getDerivativeBart0[alpha, phi, N0, n]; score = ((p1/p)[[#1 + 1]] & ) /@ data; {Plus @@ loglik, Plus @@ score}] (* ------ end of Newton-Raphson ---------- *) (* ------ golden search --------- *) likely0[x_, a_, phi_, N0_] := Module[{p, n = Max[x]}, p = pmfBartlett0[a, phi, n, N0]; Plus @@ Log[p[[x + 1]]]] goldenBartlett0[x_, N0_, phi_, opts___] := Module[{phi, lik, mom, maxIter, initTrip, tol, showIter}, initTrip = InitialTriplet /. {opts} /. Options[goldenBartlett]; If[initTrip == {1, 1, 1}, mom = bartlettMean[x, N0, phi]; initTrip = {mom/10, mom, 10*mom}]; If[!( (lik[initTrip[[1]]]>lik[initTrip[[2]]])&& (lik[initTrip[[3]]]>lik[initTrip[[2]]]) ), Print["Invalid initial triplet ..."]; Return[] ]; maxIter = MaxIterations /. {opts} /. Options[goldenBartlett]; tol = Tolerance /. {opts} /. Options[goldenBartlett]; showIter = ShowIterations /. {opts} /. Options[goldenBartlett]; lik = Function[{z}, -likely0[x, z, phi, N0]]; golden[lik, initTrip, maxIter, tol, showIter]] (* ----------- END of golden section search for MLE ------- *) (* =========== plating efficiency, added April 3, 2008 ============= *) newtonLDPlating0[data_List, e_, opts___] := Module[{m0, m1, score, info, showIter, tol, maxIter, initGuess,i, eta}, initGuess = InitialM /. {opts} /. Options[newtonLDPlating]; showIter = ShowIterations /. {opts} /. Options[newtonLDPlating]; maxIter = MaxIterations /. {opts} /. Options[newtonLDPlating]; tol = Tolerance /. {opts} /. Options[newtonLDPlating]; If[initGuess > 0, m0 = N[initGuess], If[Min[data]==0, m0 = P0Plating[data, e], m0=JonesMedianPlating[data,e]]]; If[showIter, Print[m0]]; eta=etaSeq[e, Max[data]]; For[i = 1, i <= maxIter, i++, {score, info} = LDScore$Info4Plating0[m0, e, eta, data]; m1 = m0 + score/info; If[Abs[(m1 - m0)/m0] < tol, Return[m1], If[showIter, Print[m1]]; m0 = m1; ]]; Print["no convergence!"]; ] (* February 28, 2007: to speed up *) newtonLDPlating[data_List, e_, opts___] := Module[{m0, m1, score, info, showIter, tol, maxIter, initGuess,i, eta}, initGuess = InitialM /. {opts} /. Options[newtonLDPlating]; showIter = ShowIterations /. {opts} /. Options[newtonLDPlating]; maxIter = MaxIterations /. {opts} /. Options[newtonLDPlating]; tol = Tolerance /. {opts} /. Options[newtonLDPlating]; If[initGuess > 0, m0 = N[initGuess], If[Min[data]==0, m0 = P0Plating[data, e], m0=JonesMedianPlating[data,e]]]; If[showIter, Print[m0]]; eta=etaSeq[e, Max[data]]; For[i = 1, i <= maxIter, i++, {score, info} = LDScore$Info4Plating[m0, e, eta, data]; m1 = m0 + score/info; If[Abs[(m1 - m0)/m0] < tol, Return[m1], If[showIter, Print[m1]]; m0 = m1; ]]; Print["no convergence!"]; ] (* ---- empirical mle as proposed by Lea and Coulson, extended by Angerer --- *) empiricalMLE[data_, e_, opts___] := Module[ {m0, m1, score, info, i, initGuess, maxIter, showIter, tol}, initGuess = InitialM /. {opts} /. Options[empiricalMLE]; maxIter = MaxIterations /. {opts} /. Options[empiricalMLE]; showIter = ShowIterations /. {opts} /. Options[empiricalMLE]; tol = Tolerance /. {opts} /. Options[empiricalMLE]; m0 = initGuess; If[showIter, Print[m0]]; For[i = 1, i <= maxIter, i++, score = Plus @@ (angScore[#1, m0, e] & ) /@ data; info = Plus @@ (angFisher[#1, m0, e] & ) /@ data; m1 = m0 - score/info; If[Abs[(m1 - m0)/m0] < tol, Return[m1], If[showIter, Print[m1]]; m0 = m1]; ]; Print["no convergence!"]; ] (* MathLink version for the pmf of LD with plating, 2-28-2007 *) pmfLDPlating[m_, e_, n_] := $luria`$pmfLDPlat[N[m], e, etaSeq[e, n]] pmfLDPlating0[m_, e_, n_] := Module[{k, p, q, i, j}, k=e/(1-e); p=Table[0,{n+1}]; q=etaSeq[e,n]; q = q*m*(e/(1 - e)); p[[1]] = Exp[q[[1]]]; For[i = 1, i <= n, i++, For[j = 1, j <= i, j++, p[[i + 1]] = p[[i + 1]] + j*q[[j + 1]]*p[[i - j + 1]]]; p[[i + 1]] = p[[i + 1]]/i]; p] (* ---------------- confidence interval ------------------------ *) CILDPlating0[data_List, e_, opts___] := Module[{mhat, la, qa, h, m0, m1, like, score, info, mLow, mUp, alpha, tol, initM, initLow, initUp, maxIter, showIter, i, eta}, alpha = Alpha /. {opts} /. Options[CILDPlating]; tol = Tolerance /. {opts} /. Options[CILDPlating]; maxIter = MaxIterations /. {opts} /. Options[CILDPlating]; initM = InitialM /. {opts} /. Options[CILDPlating]; initLow = InitialLowerLimit /. {opts} /. Options[CILDPlating]; initUp = InitialUpperLimit /. {opts} /. Options[CILDPlating]; showIter = ShowIterations /. {opts} /. Options[CILDPlating]; If[e <= 0 || e >= 1, Print[ "Plating efficiency \[Epsilon] out of range!"]]; eta=etaSeq[e, Max[data]]; mhat = newtonLDPlating0[data, e, InitialM -> initM, Tolerance -> tol]; {like, score} = LDLoglikely$Score4Plating0[mhat, e, eta, data]; {score, info} = LDScore$Info4Plating0[mhat, e, eta, data]; qa = Quantile[ChiSquareDistribution[1], 1 - alpha]; h = Sqrt[qa/info]; la = like - 0.5*qa; If[showIter, Print["Iterating for lower limit ..."]]; If[initLow <= 0, m0 = mhat - 0.5*h, m0 = initLow]; For[i = 1, i <= maxIter, i++, {like, score} = LDLoglikely$Score4Plating0[m0, e, eta, data]; m1 = m0 - (like - la)/score; If[showIter, Print[m0]]; If[Abs[(m1 - m0)/m0] < tol, mLow = m1; Break[], m0 = m1]]; If[showIter, Print["Iterating for upper limit ..."]]; If[initUp <= 0, m0 = mhat + 0.5*h, m0 = initUp]; For[i = 1, i <= maxIter, i++, {like, score} = LDLoglikely$Score4Plating0[m0, e, eta, data]; m1 = m0 - (like - la)/score; If[showIter, Print[m0]]; If[Abs[(m1 - m0)/m0] < tol, mUp = m1; Break[], m0 = m1]]; {mLow, mUp}] (* ------------ confidence interval, faster, March 1, 2007 ------------- *) CILDPlating[data_List, e_, opts___] := Module[{mhat, la, qa, h, m0, m1, like, score, info, mLow, mUp, alpha, tol, initM, initLow, initUp, maxIter, showIter, i, eta}, alpha = Alpha /. {opts} /. Options[CILDPlating]; tol = Tolerance /. {opts} /. Options[CILDPlating]; maxIter = MaxIterations /. {opts} /. Options[CILDPlating]; initM = InitialM /. {opts} /. Options[CILDPlating]; initLow = InitialLowerLimit /. {opts} /. Options[CILDPlating]; initUp = InitialUpperLimit /. {opts} /. Options[CILDPlating]; showIter = ShowIterations /. {opts} /. Options[CILDPlating]; If[e <= 0 || e >= 1, Print[ "Plating efficiency \[Epsilon] out of range!"]]; eta=etaSeq[e, Max[data]]; mhat = newtonLDPlating[data, e, InitialM -> initM, Tolerance -> tol]; {like, score} = LDLoglikely$Score4Plating[mhat, e, eta, data]; {score, info} = LDScore$Info4Plating[mhat, e, eta, data]; qa = Quantile[ChiSquareDistribution[1], 1 - alpha]; h = Sqrt[qa/info]; la = like - 0.5*qa; If[showIter, Print["Iterating for lower limit ..."]]; If[initLow <= 0, m0 = mhat - 0.5*h, m0 = initLow]; For[i = 1, i <= maxIter, i++, {like, score} = LDLoglikely$Score4Plating[m0, e, eta, data]; m1 = m0 - (like - la)/score; If[showIter, Print[m0]]; If[Abs[(m1 - m0)/m0] < tol, mLow = m1; Break[], m0 = m1]]; If[showIter, Print["Iterating for upper limit ..."]]; If[initUp <= 0, m0 = mhat + 0.5*h, m0 = initUp]; For[i = 1, i <= maxIter, i++, {like, score} = LDLoglikely$Score4Plating[m0, e, eta, data]; m1 = m0 - (like - la)/score; If[showIter, Print[m0]]; If[Abs[(m1 - m0)/m0] < tol, mUp = m1; Break[], m0 = m1]]; {mLow, mUp}] (* -------- Jones Median Estimator with imperfect plating efficiency ------- *) JonesMedianPlating[data_, e_] := Module[{r}, r = N[Median[data]]; If[r==0, Print["Because the sample median is zero, Jones' median method is not applicable."]; Return[] ]; (r/e - Log[2])/(Log[r/e] - Log[Log[2]])] (* ----------- Angerer Median estimator based on Lea and Coulson -------- *) LCAMedian[data_, e_, opts___] := Module[{med, m0, m1, showIter, initGuess, tol, maxIter, i}, initGuess = InitialM /. {opts} /. Options[LCAMedian]; showIter = ShowIterations /. {opts} /. Options[LCAMedian]; maxIter = MaxIterations /. {opts} /. Options[LCAMedian]; tol = Tolerance /. {opts} /. Options[LCAMedian]; med = Median[data]; If[med==0, Print["Because the sample median is zero, the method is not applicable."]; Return[] ]; m0 = initGuess; If[showIter, Print[m0]]; For[i = 1, i < maxIter, i++, m1 = m0 - (med/e/m0 - Log[m0] - 1.24)/(-med/e/m0^2 - 1/m0); If[Abs[(m1 - m0)/m0] < tol, Return[m1], If[showIter, Print[m1]]; m0 = m1]]; Print["no convergence!"]; ] (* ------------- psudo-MLE by Angerer --------------- *) mleAngerer[data_, e_, try_:-1] := Module[{mle, m, lik, guess}, If[try == -1, guess = LCAMedian[data, e], guess = try]; lik = (Exp[-0.5*(-2.02 + 11.6/(4.5 + #1/(e*m) - Log[m]))^2]/ (e*m*(4.5 + #1/(e*m) - Log[m])^2) & ) /@ data; lik = Plus @@ Log[lik]; mle = FindMinimum[-lik, {m, guess}]; mle[[2]][[1]][[2]] ] (* --------- Modified P0 method ----------------*) P0Plating[data_List, e_] := Module[{n, z}, If[e <= 0 || e >= 1, Return["Plating efficiency out of range."]]; z = Count[data, 0]; If[z == 0, Return["P0 method is not applicable."]]; n = Length[data]; (1 - e)*(Log[z/n]/e/Log[e])] (* ---- synchronized model due to Nadas and Natarajan ----- *) Nadas[data_, Nt_, N0_] := Module[{g, m}, g = Log[2, Nt/N0]; m = Mean[data]; 2*(1. - (1 - m/Nt)^(1/g))] (* ---------------- private functions ----------------------- *) (* computing the eta sequence *) etaSeq[e_, n_] := If[e <= 0.5, etaSmall[e, n], etaBig[e, n]] etaSmall[e_, n_] := Module[{odds, eta, i}, odds = e/(1 - e); eta = Table[0, {n}]; eta[[1]] = -1 - Log[e]/(1 - e); For[i = 1, i <= n - 1, i++, eta[[i + 1]] = (-odds)*eta[[i]] + 1/i/(i + 1)]; Prepend[eta, Log[e]]] etaBig[e_, n_] := Module[{revOR, eta, i}, revOR = (1 - e)/e; eta = Table[0, {n}]; eta[[n]] = ((1 - e)/n/(n + 1))*Hypergeometric2F1[1, 2, n + 2, 1 - e]; For[i = n - 1, i > 0, i--, eta[[i]] = revOR*(1/i/(i + 1) - eta[[i + 1]])]; Prepend[eta, Log[e]]] (* --------- end of the eta sequence -------------- *) (* February 28, 2007 *) seqConv[x_List, y_List]:=$luria`$pdfconv[x,y] seqConv0[x_List, y_List] := Module[{z, n = Length[x]}, z = (ListConvolve[x[[Range[#1]]], y[[Range[#1]]]] & ) /@ Range[n]; Flatten[z]] averageMutantPlating[e_,n_]:=If[e<=0.5,averageMutantPSmall[e,n], averageMutantPBig[e,n]] LDLoglikely$Score4Plating0[m_, e_, eta_List, data_List] := Module[ {k, n, p, p1, loglike, score}, k = e/(1 - e); n = Max[data]; p = pmfLDPlating0[m, e, n]; p1 = seqConv0[eta, p]*k; loglike = (Log[p[[#1 + 1]]] & ) /@ data; score = ((p1/p)[[#1 + 1]] & ) /@ data; {Plus @@ loglike, Plus @@ score}] (* acceleration by MathLink, March 1, 2007 *) LDLoglikely$Score4Plating[m_, e_, eta_List, data_List] := Module[ {k, n, p, p1, loglike, score}, k = e/(1 - e); n = Max[data]; p = pmfLDPlating[m, e, n]; p1 = seqConv[eta, p]*k; loglike = (Log[p[[#1 + 1]]] & ) /@ data; score = ((p1/p)[[#1 + 1]] & ) /@ data; {Plus @@ loglike, Plus @@ score}] LDScore$Info4Plating0[m_, e_, eta_List, data_List] := Module[{odds, n, p, p1, p2, score, info}, odds = e/(1 - e); n = Max[data]; p = pmfLDPlating0[m, e, n]; p1 = seqConv0[eta, p]*odds; p2 = seqConv0[eta, p1]*odds; score = ((p1/p)[[#1 + 1]] & ) /@ data; info = (((p1/p)^2 - p2/p)[[#1 + 1]] & ) /@ data; {Plus @@ score, Plus @@ info}] LDScore$Info4Plating[m_, e_, eta_List, data_List] := Module[{k, n, p, p1, p2, score, info}, k = e/(1 - e); n = Max[data]; p = pmfLDPlating[m, e, n]; p1 = seqConv[eta, p]*k; p2 = seqConv[eta, p1]*k; score = ((p1/p)[[#1 + 1]] & ) /@ data; info = (((p1/p)^2 - p2/p)[[#1 + 1]] & ) /@ data; {Plus @@ score, Plus @@ info}] (* ------ needed by empiricalMLE ------ *) angScore[x_, m_, e_] := (-79.741*e^3*m^3 - 7.934*e^2*m^2*x + 29.932*e*m*x^2 + x^3 + e*m*Log[m]*(19.318*e^2*m^2 - 18.432*e*m*x - x^2 + e*m*Log[m]*(-11.5*e*m - x + e*m*Log[m])))/ (m*(4.5*e*m + x - e*m*Log[m])^3) angFisher[x_, m_, e_] := (206.542*e^4*m^4 - 175.504*e^3*m^3*x - 345.152*e^2*m^2*x^2 - 16*e*m*x^3 - x^4 + e*m*Log[m]*(-231.536*e^3*m^3 + 119.792*e^2*m^2*x + 80.864*e*m*x^2 + 4*x^3 + e*m*Log[m]*(73.068*e^2*m^2 - 48.864*e*m*x - 4*x^2 + e^2*m^2*Log[m]*(-16 + Log[m]))))/(m^2*(4.5*e*m + x - e*m*Log[m])^4) (* ======== end of plating =============== *) (* ======= B0 distribution, added April 6, 2008 ======= *) pmfBartzero[A_,k_,n_]:=$luria`$pmfBartzero[A,k,n] (* Feb 11, 2008 *) (* ------------ the B0 distribution, February 15, 2008 ------------- *) (* ----------- simulating the B0 distribution, Feb 15, 2008 ----------- *) simuBartzero[A_, k_, n_, M_Integer] := Module[{prob, u, x, cumu, i, j}, prob = $luria`$pmfBartzero[A, k, M]; cumu = FoldList[Plus, First[prob], Rest[prob]]; x = Table[-1, {n}]; For[j = 1, j <= n, j++, u = Random[]; For[i = 0, i <= M, i++, If[u <= cumu[[i + 1]], x[[j]] = i; Break[]]]; ]; x] (* ------- newton for B0(A,k), Feb 15, 2008 ----------- *) (* -------- needs listConv --------- *) newtonBartzero[data_List, k_, opts___] := Module[{A0, A1, score, info, showIter, tol, maxIter, initGuess, i}, initGuess = InitialA /. {opts} /. Options[newtonBartzero]; showIter = ShowIterations /. {opts} /. Options[newtonBartzero]; maxIter = MaxIterations /. {opts} /. Options[newtonBartzero]; tol = Tolerance /. {opts} /. Options[newtonBartzero]; If[initGuess > 0, A0 = initGuess, A0 = p0Bartzero[data, k]]; If[showIter, Print[A0]]; For[i = 1, i <= maxIter, i++, {score, info} = scoreBartzero[A0, k, data]; A1 = A0 + score/info; If[Abs[(A1 - A0)/A0] < tol, Return[A1], If[showIter, Print[A1]]; A0 = A1; ] ]; Print["no convergence!"]; ] p0Bartzero[data_, k_] := Module[{p0}, If[Min[data] > 0, Return["\!\(P\_0\) method not applicable!"]]; p0 = Count[data, 0]/Length[data]; p0^(-k^(-1)) - 1] scoreBartzero[A_, k_, data_List] := Module[{n, h, p, p1, p2, score, info}, n = Max[data]; p = $luria`$pmfBartzero[A, k, n]; {p1, p2} = getBartzeroDeriv[A, k, n]; score = ((p1/p)[[#1 + 1]] & ) /@ data; info = (((p1/p)^2 - p2/p)[[#1 + 1]] & ) /@ data; {Plus @@ score, Plus @@ info}] getBartzeroDeriv[A_, k_, n_] := Module[{p1, p2, h, h2, j}, h = Prepend[Table[1.0/(j*(j + 1)), {j, 1, n}], -1]; p1 = k*$luria`$pdfconv[$luria`$pmfBartzero[A, k + 1, n], h]; h2 = $luria`$pdfconv[h, h]; p2 = k*(k + 1)*$luria`$pdfconv[$luria`$pmfBartzero[A, k + 2, n], h2]; {p1, p2}] (* ------- CI for B0(A,k), Feb 22, 2008 ----------- *) CIBartzero[data_List, k_, opts___] := Module[{Ahat, la, qa, h, A0, A1, like, score, info, ALow, AUp, alpha, tol, initA, initLow, initUp, maxIter, showIter, i}, alpha = Alpha /. {opts} /. Options[CIBartzero]; tol = Tolerance /. {opts} /. Options[CIBartzero]; showIter = ShowIterations /. {opts} /. Options[CIBartzero]; maxIter = MaxIterations /. {opts} /. Options[CIBartzero]; initA = InitialA /. {opts} /. Options[CIBartzero]; initLow = InitialLowerLimit /. {opts} /. Options[CIBartzero]; initUp = InitialUpperLimit /. {opts} /. Options[CIBartzero]; If[showIter, Print["Iterating for mle of A..."]]; Ahat = newtonBartzero[data, k, InitialA -> initA, ShowIterations->showIter]; {like, score} = loglikScoreBartzero[Ahat, k, data]; {score, info} = scoreBartzero[Ahat, k, data]; qa = Quantile[ChiSquareDistribution[1], 1 - alpha]; h = Sqrt[qa/info]; la = like - 0.5*qa; If[showIter, Print["Iterating for lower limit ..."]]; If[initLow <= 0, A0 = Ahat - 0.5*h, A0 = initLow]; For[i = 1, i <= maxIter, i++, {like, score} = loglikScoreBartzero[A0, k, data]; A1 = A0 - (like - la)/score; If[showIter, Print[A0]]; If[Abs[(A1 - A0)/A0] < tol, ALow = A1; Break[], A0 = A1]]; If[showIter, Print["Iterating for upper limit ..."]]; If[initUp <= 0, A0 = Ahat + 0.5*h, A0 = initUp]; For[i = 1, i <= maxIter, i++, {like, score} = loglikScoreBartzero[A0, k, data]; A1 = A0 - (like - la)/score; If[showIter, Print[A0]]; If[Abs[(A1 - A0)/A0] < tol, AUp = A1; Break[], A0 = A1]]; {ALow, AUp}] (* needed only by CIBartzero, 2-22-08 *) loglikScoreBartzero[A_, k_, data_List] := Module[{n, h, p, p1, p2, score, loglik}, n = Max[data]; p = $luria`$pmfBartzero[A, k, n]; loglik = (Log[p[[#1 + 1]]] & ) /@ data; {p1, p2} = getBartzeroDeriv[A, k, n]; score = ((p1/p)[[#1 + 1]] & ) /@ data; {Plus @@ loglik, Plus @@ score}] (* ---------- end of B0 distribution ------------- *) End[] EndPackage[]