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

[John Denker <jsd@av8n.com>]

On 9/22/21 8:14 AM, Paul Nord wrote:

> Traditionally we have them fit a double exponential function to this data.
> But this is fraught with problems and subtle choices in the analysis which
> will change the result.  And it's not how these half lives are really
> measured.  Though, obviously in a precision measurement you'd want to have
> continuous data collection.

I assume the number of counts in any given interval is not large.
Otherwise you wouldn't be asking the question. Generic least-squares
fitting would do the job. Crude but effective.

Core problem #1 is that "least squares" routines assume the statistics
are Gaussian. The logarithm of a Gaussian is a quadratic. Adding the
quadratics is isomorphic to multiplying the probabilities to find the
joint probability. Hence least squares.

That's a problem because when the counts are small, a Poisson distribution
doesn't resemble a Gaussian. Large counts yes, small counts no.

Core problem #2 is the time dependence. You can find a squillion youtube
videos on how to find *the* rate of a Poisson process when the rate isn't
changing.

===========

As a warmup exercise, let's temporarily switch from analysis of data to
synthesis of data. That, suppose we wanted to simulate the experiment.
Given the four parameters N1(initial), N2(initial), τ1, and τ2, you
can figure out the Poisson rate λ at any given time. Then for any given
interval near that time, you can synthesize an appropriate number of
events for that interval.

So that's the model. Now the task is to adjust the parameters of the
model to fit the real observed data.

Given a set of parameters, then for each observation it is straightforward
to calculate the probability that the model, using those parameters, would
have emitted the observed event-count.  Take the log of the probability.
Sum over all observations. That's the objective function. Search the
parameter space to optimize the objective function.

That's it. That's a decent way to formalize the problem.

Now it's "just" a four-dimensional search problem. Pain in the neck.
If necessary you can stabilize it by first doing a 2D fit to estimate
the fast exponential, then a 2D fit for the slow one, then a joint 4D
fit.

I realize this does not answer the question, insofar as it is not
particularly Bayesian. It is still maximum_a_priori (boo!) i.e. the
conditional probability of the data given the parameters ... when
what you really want is maximum_a_posteriori i.e. the probability
of the parameters given the data.

You could fix that given a Bayesian prior on the parameters, but at
the moment I have no clue what the prior should be.

In the curve fitting business the priors often take the form of
"regularizers" but it's not obvious what regularizers to use. I
suppose you could constrain the search to keep the lifetimes within
(say) a factor of 2 of the book value......

Given a decent amount of /total/ data the priors shouldn't matter
too much. The total number of events can be large even if the number
of events in any particular observation is not.

If the total number of events is tiny then a reliable analysis may
be simply impossible.