| 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] |
Oops! My solution is no good because my log-likelihood function uses, No,
the number of radioactive atoms in the sample at time t = 0. For the
correct Baysian solution, I will have to find the joint posterior
probability for the half-life and for, No, and then marginalize out No,
i.e. integrate over No to get a univariate posterior probability for tau
alone....
stay tuned. :-)
On Tue, Sep 21, 2021 at 11:30 PM Francois Primeau <fprimeau@uci.edu>
wrote:
I'm assuming the problem is to estimate the half-life of a sample given
count data taken at different times. To be concrete I'm going to assume
that you have counts, c1, c2, and c3 taken at times t1, t2, and t3. Also,
I'm going to assume that the counting was done for 10 minutes, so
dt1=dt2=dt3 = 10*60 seconds. I'm going to assume that the error in the
data is negligible and that only the Poisson statistics matter. The Matlab
code and the comments in it explain the solution. The code starts off by
generating synthetic data and then using the synthetic data to estimate the
half-life.
The Bayesian solution (given in Matlab code) is the following:
tau_true = 52*60*60; % half-life of neutron activated copper in seconds
No = 1e6; % 1 million radioactive atoms at t0 = 0
% all times are in seconds
t1 = 0; dt1 = 10*60; % count for 10 minutes at t1
t2 = 24*60^2; dt2 = 10*60; % count for 10 minutes at t2
t3 = 48*60^2; dt3 = 10*60; % count for 10 minutes at t3
% generate synthetic data
c1 = poissrnd( log(2)*No*exp(-log(2)*t1/tau_true)*dt1/tau_true )
c2 = poissrnd( log(2)*No*exp(-log(2)*t2/tau_true)*dt2/tau_true )
c3 = poissrnd( log(2)*No*exp(-log(2)*t3/tau_true)*dt3/tau_true )
% Now pretend we don't know tau and want to use Bayes' theorem to
estimate it
% Poisson probability function: prob(c) = mu^c*exp(-mu)/c!
% log likelihood function for one of the c's (+ an irrelevant constant
that doesn't depend on tau )
ll = @(tau,c,t,dt) c * log( log(2)*No*exp(-log(2)*t./tau)*dt./tau ) -
log(2)*No*exp(-log(2)*t./tau)*dt./tau;
% The joint probability factors because the probabilities for the counts
*conditioned on tau* are independent
% so the log-likelihood for the full data set is given by the sum of the
individual log-likelihoods
% log-likelihood for the full data set:
LL = @(tau) ll(tau,c1,t1,dt1) + ll(tau,c2,t2,dt2) + ll(tau,c3,t3,dt3);
% I'm playing dumb and assuming the only thing I know is that tau is a
% positive number. In practice, I could use a more informative prior,
% but if I collect even a moderate amount of good data, then the prior
won't
% matter much.
% prior: prob(log(tau)) = constant ==> prob(tau) propto 1/tau
lp = @(tau) 1./tau; % log prior (not normalized!)
% log posterior from Bayes' Theorem (not normalized so there is a missing
constant)
LP = @(tau) LL(tau)+lp(tau);
% now make a plot of the posterior
tau = linspace(1.5e5,2.5e5,100); % s
P = exp(LP(tau)-max(LP(tau))); % prevent overflow
Z = trapz(tau,P); % normalization
P = P/Z; %normalized posterior
plot(tau/60,P*60); grid on
xlabel('half life (minutes)')
ylabel('probability density (1/minutes)')
% the plot shows a nice bell shaped posterior
%so the mean and variance will provide a good summary of the posterior
% compute the mean of the posterior
tau_hat = trapz(tau,P.*tau);
% compute the variance of the posterior
sig2 = trapz(tau,P.*(tau-tau_hat).^2);
% print the result
fprintf('tau = %f +/- %f minutes \n',tau_hat/60,sqrt(sig2)/60)
On Tue, Sep 21, 2021 at 10:29 AM Paul Nord <Paul.Nord@valpo.edu> wrote:
That title doesn't google well. It seems like there should be a good
reference for this. The best I've got is David MacKay's code here:
http://www.inference.org.uk/mackay/itprnn06/slides/10/mgp00010.html
Translated into python below.
Two questions:
1. If we record the number of counts observed with a geiger tube at
discrete periods within the decay, is this approach still valid? Say
that
I've got a sample with a 52 hour half life. I come back about once a
day,
fire up the good old geiger tube and measure the activity for 10 minutes.
Can I just use that number of counts as the power of a probability
function
to multiply through here?
2. Will this give me a good result if I'm extracting multiple decay
periods
from the same data?
The half lives of neutron activated copper is the experiment, actually.
Paul
import matplotlib.pyplot as plt
import numpy as np
import math
a = 1.0
b = 20.0
def p(x,l):
val = math.exp( -x/l ) / Z(l)/l
return val
def Z(l):
val = math.exp( -a/l) - math.exp(-b/l)
return val
def like(l):
return p(3,l) * p(5,l) * p(12,l)
t1 = np.arange(0.01,1000.0,.1)
fig,axs = plt.subplots(2)
axs[0].semilogx()
axs[1].semilogx()
axs[0].plot(t1,[p(3,x) for x in t1])
axs[0].plot(t1,[p(5,x) for x in t1])
axs[0].plot(t1,[p(12,x) for x in t1])
axs[1].plot(t1,[like(x) for x in t1])
plt.show()
_______________________________________________
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
--
Professor of Earth System Science
Rm 3224 Croul Hall
University of California, Irvine