{"id":674,"date":"2021-09-05T18:13:38","date_gmt":"2021-09-05T18:13:38","guid":{"rendered":"https:\/\/www.allendowney.com\/blog\/?p=674"},"modified":"2021-11-19T03:55:17","modified_gmt":"2021-11-19T03:55:17","slug":"emitter-detector-redux","status":"publish","type":"post","link":"https:\/\/www.allendowney.com\/blog\/2021\/09\/05\/emitter-detector-redux\/","title":{"rendered":"Emitter Detector Redux"},"content":{"rendered":"\n

In the first edition of Think Bayes, I presented what I called the Geiger counter problem<\/a>, which is based on an example in Jaynes, Probability Theory<\/em>. But I was not satisfied with my solution or the way I explained it, so I cut it from the second edition.<\/p>\n\n\n\n

I am re-reading Jaynes now, following the excellent series of videos by Aubrey Clayton<\/a>, and this problem came back to haunt me. On my second attempt, I have a solution that is much clearer, and I think I can explain it better.<\/p>\n\n\n\n

I’ll outline the solution here, but for all of the details, you can read the bonus chapter<\/a>, or click here to run the notebook on Colab<\/a>.<\/p>\n\n\n\n

The Emitter-Detector Problem<\/h1>\n\n\n\n

Here\u2019s the example from Jaynes, page 168:<\/p>\n\n\n\n

We have a radioactive source \u2026 which is emitting particles of some sort \u2026 There is a rate p<\/em>, in particles per second, at which a radioactive nucleus sends particles through our counter; and each particle passing through produces counts at the rate \u03b8<\/em>. From measuring the number {c1 , c2 , \u2026} of counts in different seconds, what can we say about the numbers {n1 , n2 , \u2026} actually passing through the counter in each second, and what can we say about the strength of the source?<\/p><\/blockquote>\n\n\n\n

As a model of the source, Jaynes suggests we imagine \u201cN<\/em> nuclei, each of which has independently the probability r<\/em> of sending a particle through our counter in any one second\u201d. If N<\/em> is large and r<\/em> is small, the number of particles emitted in a given second is well modeled by a Poisson distribution with parameter s<\/em>=Nr<\/em>, where s<\/em> is the strength of the source.<\/p>\n\n\n\n

As a model of the sensor, we\u2019ll assume that \u201ceach particle passing through the counter has independently the probability \u03d5<\/em> of making a count\u201d. So if we know the actual number of particles, n<\/em>, and the efficiency of the sensor, \u03d5<\/em>, the distribution of the count is Binomial(n<\/em>,\u03d5<\/em>).<\/p>\n\n\n\n

With that, we are ready to solve the problem. Following Jaynes, I\u2019ll start with a uniform prior for s<\/em>, over a range of values wide enough to cover the region where the likelihood of the data is non-negligible. To represent distributions, I’ll use the Pmf class from empiricaldist.<\/p>\n\n\n\n

ss = np.linspace(0, 350, 101)\nprior_s = Pmf(1, ss)<\/pre>\n\n\n\n
For each value of s<\/em>, the distribution of n<\/em> is Poisson, so we can form the joint prior of s<\/em> and n<\/em> using the poisson<\/code> function from SciPy. The following function creates a Pandas DataFrame<\/code> that represents the joint prior.<\/p>\n\n\n\n
def make_joint(prior_s, ns):\n    ss = prior_s.qs\n    S, N = np.meshgrid(ss, ns)\n    ps = poisson(S).pmf(N) * prior_s.ps\n    joint = pd.DataFrame(ps, index=ns, columns=ss)\n    joint.index.name = 'n'\n    joint.columns.name = 's'\n    return joint<\/pre>\n\n\n\n
The result is a DataFrame with one row for each value of n<\/em> and one column for each value of s<\/em>.<\/p>\n\n\n\n
To update the prior, we need to compute the likelihood of the data for each pair of parameters. However, in this problem the likelihood of a given count depends only on n<\/em>, regardless of s<\/em>, so we only have to compute it once for each value of n<\/em>. Then we multiply each column in the prior by this array of likelihoods. The following function encapsulates this computation, normalizes the result, and returns the posterior distribution.<\/p>\n\n\n\n
def update(joint, phi, c):\n    ns = joint.index\n    likelihood = binom(ns, phi).pmf(c)\n    posterior = joint.multiply(likelihood, axis=0)\n    normalize(posterior)\n    return posterior<\/pre>\n\n\n\n
First update<\/h2>\n\n\n\n
Let\u2019s test the update function with the first example, on page 178 of Probability Theory<\/em>:<\/p>\n\n\n\n
During the first second, c1 = 10<\/code> counts are registered. What can [we] say about the number n1<\/code> of particles?<\/p><\/blockquote>\n\n\n\n
Here\u2019s the update:<\/p>\n\n\n\n
c1 = 10\nphi = 0.1\nposterior = update(joint, phi, c1)<\/pre>\n\n\n\n
The following figures show the posterior marginal distributions of s<\/em> and n<\/em>.<\/p>\n\n\n\n
posterior_s = marginal(posterior, 0)<\/pre>\n\n\n\n
\"_images\/radiation_24_1.png\"\/<\/figure>\n\n\n\n
posterior_n = marginal(posterior, 1)<\/pre>\n\n\n\n
\"_images\/radiation_25_1.png\"\/<\/figure>\n\n\n\n
The posterior mean of n<\/em> is close to 109, which is consistent with Equation 6.116. The MAP is 99, which is one less than the analytic result in Equation 6.113, which is 100. It looks like the posterior probabilities for 99 and 100 are the same, but the floating-point results differ slightly.<\/p>\n\n\n\n
Jeffreys prior<\/h2>\n\n\n\n
Instead of a uniform prior for s<\/em>, we can use a Jeffreys prior, in which the prior probability for each value of s<\/em> is proportional to 1\/s<\/em>. This has the advantage of \u201cinvariance under certain changes of parameters\u201d, which is \u201cthe only correct way to express complete ignorance of a scale parameter.\u201d However, Jaynes suggests that it is not clear \u201cwhether s<\/em> can properly be regarded as a scale parameter in this problem.\u201d Nevertheless, he suggests we try it and see what happens. Here\u2019s the Jeffreys prior for s<\/em>.<\/p>\n\n\n\n
prior_jeff = Pmf(1\/ss[1:], ss[1:])<\/pre>\n\n\n\n
We can use it to compute the joint prior of s<\/em> and n<\/em>, and update it with c1<\/code>.<\/p>\n\n\n\n
joint_jeff = make_joint(prior_jeff, ns)\nposterior_jeff = update(joint_jeff, phi, c1)<\/pre>\n\n\n\n
Here\u2019s the marginal posterior distribution of n<\/em>:<\/p>\n\n\n\n
posterior_n = marginal(posterior_jeff, 1)<\/pre>\n\n\n\n
\"_images\/radiation_35_1.png\"\/<\/figure>\n\n\n\n
The posterior mean is close to 100 and the MAP is 91; both are consistent with the results in Equation 6.122.<\/p>\n\n\n\n
Robot A<\/h2>\n\n\n\n
Now we get to what I think is the most interesting part of this example, which is to take into account a second observation under two models of the scenario:<\/p>\n\n\n\n
Two robots, [A and B], have different prior information about the source of the particles. The source is hidden in another room which A and B are not allowed to enter. A has no knowledge at all about the source of particles; for all [it] knows, \u2026 the other room might be full of little [people] who run back and forth, holding first one radioactive source, then another, up to the exit window. B has one additional qualitative fact: [it] knows that the source is a radioactive sample of long lifetime, in a fixed position.<\/p><\/blockquote>\n\n\n\n
In other words, B has reason to believe that the source strength s<\/em> is constant from one interval to the next, while A admits the possibility that s<\/em> is different for each interval. The following figure, from Jaynes, represents these models graphically.<\/p>\n\n\n\n
\"\"<\/figure>\n\n\n\n
For A, the \u201cdifferent intervals are logically independent\u201d, so the update with c2 = 16<\/code> starts with the same prior.<\/p>\n\n\n\n
c2 = 16\nposterior2 = update(joint, phi, c2)<\/pre>\n\n\n\n
Here\u2019s the posterior marginal distribution of n2<\/code>. <\/p>\n\n\n\n
\"_images\/radiation_41_1.png\"\/<\/figure>\n\n\n\n
The posterior mean is close to 169, which is consistent with the result in Equation 6.124. The MAP is 160, which is consistent with Equation 6.123.<\/p>\n\n\n\n
Robot B<\/h2>\n\n\n\n
For B, the \u201clogical situation\u201d is different. If we consider s<\/em> to be constant, we can \u2013 and should! \u2013 take the information from the first update into account when we perform the second update. We can do that by using the posterior distribution of s<\/em> from the first update to form the joint prior for the second update, like this:<\/p>\n\n\n\n
joint = make_joint(posterior_s, ns)\nposterior = update(joint, phi, c2)\nposterior_n = marginal(posterior, 1)<\/pre>\n\n\n\n
\"_images\/radiation_45_1.png\"\/<\/figure>\n\n\n\n
The posterior mean of n<\/em> is close to 137.5, which is consistent with Equation 6.134. The MAP is 132, which is one less than the analytic result, 133. But again, there are two values with the same probability except for floating-point errors.<\/p>\n\n\n\n
Under B\u2019s model, the data from the first interval updates our belief about s<\/em>, which influences what we believe about n2<\/code>.<\/p>\n\n\n\n
Going the other way<\/h2>\n\n\n\n
That might not seem surprising, but there is an additional point Jaynes makes with this example, which is that it also works the other way around: Having seen c2<\/code>, we have more information about s<\/em>, which means we can \u2013 and should! \u2013 go back and reconsider what we concluded about n1<\/code>.<\/p>\n\n\n\n
We can do that by imagining we did the experiments in the opposite order, so<\/p>\n\n\n\n
  1. We\u2019ll start again with a joint prior based on a uniform distribution for s<\/em><\/li>
  2. Update it based on c2<\/code>,<\/li>
  3. Use the posterior distribution of s<\/em> to form a new joint prior,<\/li>
  4. Update it based on c1<\/code>, and<\/li>
  5. Extract the marginal posterior for n1<\/code>.<\/li><\/ol>\n\n\n\n
    joint = make_joint(prior_s, ns)\nposterior = update(joint, phi, c2)\nposterior_s = marginal(posterior, 0)\n\njoint = make_joint(posterior_s, ns)\nposterior = update(joint, phi, c1)\nposterior_n2 = marginal(posterior, 1)<\/pre>\n\n\n\n
    The posterior mean is close to 131.5, which is consistent with Equation 6.133. And the MAP is 126, which is one less than the result in Equation 6.132, again due to floating-point error.<\/p>\n\n\n\n
    Here\u2019s what the new distribution of n1<\/code> looks like compared to the original, which was based on c1<\/code> only.<\/p>\n\n\n\n
    \"_images\/radiation_58_0.png\"\/<\/figure>\n\n\n\n
    With the additional information from c2<\/code>:<\/p>\n\n\n\n