Re: [Phys-L] Bayesian Inference in Half-Life measurement

[Brian Whatcott <betwys1@sbcglobal.net>]

Francois,
Copper naturally occurs at isotopic ratios around  2.25 : 1Moreover, the half lives in question are dramatically different;  12.7 hours and 5.12 minutes (Paul's numbers)
One hour of observation would reduce the activity of one of the species a thousand fold.If observation were continued for the succeeding four days, such activity would be reduced by a factor greater than 1 x 10 ^-27 next day, i.e negligible.   It seems then, that the old-school method of computing a decay constant for the slow speciesfrom fours days' data and using it to correct the first hour's observations for the slow species, and for background radiation would serve well.
What am I missing?
Brian   
   - Francois Primeau via Phys-l <phys-l@mail.phys-l.org>To:phys-l@phys-l.orgCc:Francois PrimeauWed, Sep 22 at 11:02 AMFitting multiple exponential decays for the case with unknown decay rates and unknown initial number of atoms is a notoriously tricky numerical problem because it’s easy to get good fits with very different parameter values. (Acton’s classic book “Numerical Methods that Usually Work” discusses the issue.)   
   
A Bayesian approach can help in this case if you happen to have additional information on top of the counts.  If you do, you can  build this information into the prior and maybe rule out lots of the posterior parameter space.   
   
My Matlab example was for the case of a single exponential decay.  If the experimental data includes two decaying isotopes, my code can be extended to produce the posterior probability for 4 parameters.  The numerical integrations over a 4d mesh of points will get quite a bit more costly in terms of computer time, but it should still take less than a few minutes on a modern laptop.  I’m pretty sure though that the posterior probability will be quite broad and produce large error bars for the estimated decay constants unless extra info is built into the posterior. But that’s not be a problem with the estimation method. It’s a feature of the problem.   
   
Sent from my iPhone   
   
> On Sep 22, 2021, at 8:29 AM, Brian Whatcott <betwys1@sbcglobal.net> wrote:   
>   
>  Francois' Take 2 is not playing happily with the other children at present. I may have erred in the transcription.   
>     On Wednesday, September 22, 2021, 10:00:15 AM CDT, Brian Whatcott <betwys1@sbcglobal.net> wrote:    
>   
>   Like Francois, I played with a MATLAB code for exponential decay. Galloping off rapidly in all directions, 
   - I generated two exponential decay series with noise, to address the issue of extracting data for two decay species,
   -  like this.....[spoiler: extracting multiple decay parameters from observations is an on going topic in the literature]
   -    

   - >> >> x = (0:0.2:5)';y = 2*exp(-0.2*x) + 0.01*randn(size(x));y2 = 3*exp(-0.3*x) + 0.01*randn(size(x));   
> f = fit(x,y,'exp1')plot(f,x,y)   
> f =   
>       General model Exp1:    f(x) = a*exp(b*x)    Coefficients (with 95% confidence bounds):      

   -    a =        2.06  (1.941, 2.179)      b =    -0.1894  (-0.2159, -0.1629)   
>   
> > > f2 = fit(x,y2,'exp2')plot(f2,x,y2)   
> f2 =   
>       General model Exp2:    f2(x) = a*exp(b*x) + c*exp(d*x)    Coefficients (with 95% confidence bounds):      a =      -0.197  (-7.766e+04, 7.766e+04)      b =    -0.2926  (-1567, 1566)      c =      3.203  (-7.766e+04, 7.767e+04)      d =    -0.3006  (-98.32, 97.72)>>   
> [Conclusion: stop wasting time!]
   -    
> Next I plotted Francois' 'nice-shaped posterior' (forgive the mental image) whichI show here:https://imgur.com/gallery/TpP4ZdY   
>   
> ...to confirm his code for a noiseless error-free data set.   
>