Chronology Current Month Current Thread Current Date
[Year List] [Month List (current year)] [Date Index] [Thread Index] [Thread Prev] [Thread Next] [Date Prev] [Date Next]

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

@Brian: I think your intuition is correct that if the decay rates are
sufficiently different, the "old school" method would work fine, provided
you collected data multiple times early on to resolve the decay of the
first isotope and then more data that are spaced further apart to resolve
the decay of the longer-lived isotope.

Now I'm curious to try it with simulated data and background noise to see
how well the individual decay rates can be estimated as the number of
counts gets small. The problem is that the posterior probability will be
over a 5-dimensional space, and doing the simple (and transparent)
discretization of the posterior over a regular mesh will become quite
expensive. So I might need to be smarter about how I marginalize out the
nuisance parameters, i.e., how I integrate over N1o, N2o, and B in the
joint posterior for (tau1, tau2, N1o, N2o, and B) to get a joint posterior
for tau1 and tau2 alone. It makes the problem less fun.


On Wed, Sep 22, 2021 at 11:29 AM Brian Whatcott <betwys1@sbcglobal.net>
wrote:

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.



_______________________________________________
Forum for Physics Educators
Phys-l@mail.phys-l.org
https://www.phys-l.org/mailman/listinfo/phys-l



--
Professor of Earth System Science
Rm 3224 Croul Hall
University of California, Irvine