1 Set ⃗ n$_{o}$ to the array of observed event numbers. 2 For i = 1 , . . . , I : 3 Generate ( σ$_{i}$, ⃗ ϵ$_{i}$, ⃗ µ$_{i}$ ) ∼ p ( ⃗ n$_{o}$ | σ, ⃗ ϵ, ⃗ µ ) π ( ⃗ ϵ, ⃗ µ ). 4 For j = 1 , . . . , J : 5 Generate ⃗ n$_{j}$ ∼ p ( ⃗ n | σ$_{i}$, ⃗ ϵ$_{i}$, ⃗ µ$_{i}$ ). 6 Calculate d $^{2}$[ − ln p ( ⃗ n$_{j}$ | σ$_{i}$, ⃗ ϵ$_{i}$, ⃗ µ$_{i}$ )] /dσ 2 i by numerical differentiation. 7 Average the J values of d $^{2}$[ − ln p ( ⃗ n | σ$_{i}$, ⃗ ϵ$_{i}$, ⃗ µ$_{i}$ )] /dσ 2 i obtained at line 6, and take the square root. This yields a numerical approximation to the conditional Jeffreys’ prior π$_{J}$ ( σ$_{i}$ | ⃗ ϵ$_{i}$, ⃗ µ$_{i}$ ). 8 Histogram the σ$_{i}$ values generated at line 3, weighting them by π$_{J}$ ( σ$_{i}$ | ⃗ ϵ$_{i}$, ⃗ µ$_{i}$ ) /p ( ⃗ n$_{o}$ | σ$_{i}$, ⃗ ϵ$_{i}$, ⃗ µ$_{i}$ ). This yields π$_{R}$$_{1}$ ( σ ), the σ -marginal prior. 9 Histogram the σ$_{i}$ values generated at line 3, weighting them by π$_{J}$ ( σ$_{i}$ | ⃗ ϵ$_{i}$, ⃗ µ$_{i}$ ). This yields π$_{R}$$_{1}$ ( σ | ⃗ n$_{o}$ ), the σ -marginal posterior. Although not required for the calculation of the reference posterior, an approximation to the reference prior is provided at line 8. By construction this approximation is only reliable for σ values in the bulk of the flat-prior posterior. The generation step at line 3 is done via a Markov chain Monte Carlo procedure [19]. The particular choice of sampling distribution for the generated ( σ, ⃗ ϵ, ⃗ µ ) triplets is motivated by the desire to obtain weights with reasonably small variance at steps 8 and 9. However, the flat-prior posterior p ( ⃗ n$_{0}$ | σ, ⃗ ϵ, ⃗ µ ) π ( ⃗ ϵ, ⃗ µ ) is not always proper with respect to ( σ, ⃗ ϵ, ⃗ µ ). When M = 1 (single-count model), it is improper if x ≤ 1 / 2. Propriety can then be restored by multiplying the flat-prior posterior by ϵ and correspondingly adjusting the weights at steps 8 and 9. Another feature of the above algorithm is that it does not implement the compact subset normalization. In the cases that we examined, this procedure made no difference, but this may not be true for more general problems than those our code seeks to solve. Unfortunately the current lack of guidelines in the choice of compact sets limits our ability to address this issue in the code. The algorithm for Method 2 has a simpler structure, since all it does is apply Jeffreys’ rule to a marginalized likelihood p ( ⃗ n$_{o}$ | σ ) provided by the user. The calculation does not require random sampling of the parameters and is done at fixed σ values. For a given σ , the reference prior π$_{R}$$_{2}$ ( σ ) is obtained by Monte Carlo averaging, over an ensemble of vectors 17