# Bayesian non-parametric modeling by mixture of Kibble-Pólya tree to detect Low-Activity uranium contamination

Hanan Arahmane, Jonathan Dumazert, Eric Barat, Thomas Dautremer, Frédérick Carrel, Nicolas Dufour, Maugan Michel

► To cite this version:

Hanan Arahmane, Jonathan Dumazert, Eric Barat, Thomas Dautremer, Frédérick Carrel, et al. Bayesian non-parametric modeling by mixture of Kibble-Pólya tree to detect Low-Activity uranium contamination. Results in Physics, 2024, 64, pp.107874. 10.1016/j.rinp.2024.107874. cea-05449063

HAL Id: cea-05449063

https://cea.hal.science/cea-05449063v1

Submitted on 8 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Bayesian non-parametric modeling by mixture of Kibble-Pólya tree to detect Low-Activity uranium contamination

Hanan Arahmane<sup>a, *</sup>, Jonathan Dumazert<sup>b</sup>, Eric Barat<sup>a</sup>, Thomas Dautremer<sup>a</sup>, Frédérick Carrel<sup>a</sup>, Nicolas Dufour<sup>b</sup>, Maugan Michel<sup>a</sup>

<figure><img src="image_2.png" /></figure>

Check for
updates

## ARTICLE INFO

**Keywords:**

Finite Pólya tree model

Non-parametric Bayesian modeling

Gamma rays

Low-level radioactivity

Spectrometry

## ABSTRACT

Accurate detection of low-level radioactivity is critical for decommissioning projects in nuclear facilities, particularly in the design of radiation monitoring systems with a low false alarm rate. Utilizing a non-parametric Bayesian continuous probability distribution enables reliable mapping of potential contamination. Our method introduces a statistical test based on a Pólya tree prior, applied to radiation detection. The detection efficiency of this proposed Bayesian test is quantified using receiver-operating characteristic (ROC) curves and compared to a Bayesian test based on the Kibble bivariate gamma distribution developed for the same purpose. The results demonstrate that the new Bayesian test generally outperforms the previous method in terms of detection performance under very low signal-to-noise ratios, with improvements ranging from 3% to 28% against both stationary and non-stationary radiological backgrounds, respectively. This superiority is further reaffirmed through comparisons with alternative Bayesian and frequentist hypothesis tests, with gains estimated at 52% and 4%, respectively.

## Introduction

Bayesian Inference (BI) decision-making has recently become a revolutionary method in radiation detection, especially vital for decommissioning and dismantling nuclear facilities, where achieving very low detection limits is essential. BI provides a reliable and accurate framework for identifying the most effective strategies while meeting operational standards. This is critically important in nuclear engineering from both safety and risk management perspectives. Inaccurate detection can lead to over-classification, significantly raising disposal costs, or under-classification, compromising the safety of nearby residents. The growing number of decommissioning projects, each facing significant operational challenges, highlights the urgent need for innovative BI-based methods. Therefore, technical solutions must strictly adhere to rigorous requirements, particularly concerning detection limits.

Numerous significant contributions utilizing BI have been proposed to tackle specific challenges in various contexts. In one instance, a Bayesian framework was adopted to assess the activity of anthropogenic radionuclides in soils, as part of evaluating the environmental impact of decommissioning activities at the Garigliano Nuclear Power Plant [1].

This study focused on Cs-137 and U-236 activities in soil samples, analyzing their statistical distributions over space and time. The authors highlighted the Bayesian approach's superiority over frequentist methods, emphasizing its enhanced performance with small sample sizes and non-normal data distributions, its straightforward group comparisons, its flexibility in incorporating multiple prior distributions, and its ability to handle non-detectable values. Another application of BI involved estimating the remote depth information of localized low-level radioactive contaminants relevant to nuclear decommissioning [2]. Experimental results demonstrated the method's capability to detect both individual and multiple buried radioisotopes, providing accurate estimations of depth and activity with a 95% confidence interval. Additionally, a Bayesian network (BN) based on copulas was employed for processing system risk assessment by studying non-linear relationships among process components [3]. This approach, relying on conditional probabilities and joint probability densities, outperformed traditional BNs in a practical case study. In the field of industrial metrology, a modified Bayes' theorem was proposed for reliable conformity assessment [4]. Using Random Fuzzy variables to represent measurement results and metrological *a priori* knowledge, the authors

<footnote><sup>a</sup> Université Paris-Saclay, CEA, List, F-91120 Palaiseau, France</footnote>

<footnote><sup>b</sup> CEA-DAM, DIF, F-91297 Arpajon, France</footnote>

<footnote>* Corresponding author.
E-mail address: hanan.arahmane@gmail.com (H. Arahmane).</footnote>
demonstrated that their approach significantly reduced the risk of incorrect decisions compared to the standard Bayes' theorem, especially in cases of measuring instrument deviations. Shifting focus to gamma-ray spectra analysis, BI was applied to estimate uncertainty in isotope activity ratios, providing a more precise quantification of isotopes [5]. The method's effectiveness was evaluated using complex spectra obtained from experiments and Monte Carlo N-Particle Transport Code 6 (MCNP6) simulations, demonstrating accurate determination and activity ratio estimation in low-resolution spectra. Lastly, a hierarchical Bayesian model was employed for gamma radiative source localization and background radiation estimation [6]. The model exhibited acceptable accuracy without the need for discriminating background radiation estimation at each sensor location, and computational costs were minimized through the use of a simple photon transport model. In a different study [7], Bayesian inversion of borehole gamma spectrometry data was used to estimate soil radioactivity-depth profiles. The authors employed a Bayesian inversion method with a Gaussian process prior to determine the activity concentration of Cs-137 and its uncertainty in soil, using in-situ borehole gamma spectrometry datasets. Results highlighted soil density, activity distribution, error in detector position, and borehole diameter as crucial parameters for activity estimation. Moreover, the calculated Cs-137 activity concentrations agreed well with those from soil samples, confirming Bayesian inversion's utility for comprehensive probabilistic site characterization with extensive in-situ gamma spectrometry datasets.

In this study, our focus is on introducing a novel method based on Bayesian decision-making within the context of nuclear facility decommissioning and dismantling. Specifically, we aim to address the detection of surface activity of about 500 Bq.0.1 m⁻² uranium contamination on large-scale concrete surfaces. We utilize gamma-ray spectrometry, with a specific focus on identifying a potential radiological signature associated with uranium isotope, namely, uranium 235 and 238 (U-235 and U-238). High-purity germanium detector (HPGe) count measurements form the basis of our gamma-ray spectrometry approach [8]. However, the expected low count rate linked to this signature poses a detection challenge, with a signal-to-noise ratio (SNR) as low as 10⁻¹ within a 1 s acquisition time [9]. The main constraint arises from the presence of radiological background, a significant and fluctuating spectral contribution dependent on the scanned location. Variations in both intensity and spectral shape of the radiological background are expected across different spots, complicating the interpretation of count spectra. Furthermore, practical considerations necessitate that the duration of each measurement remains reasonable, facilitating the inspection of extensive areas that may encompass several hundreds or thousands of square meters.

In our prior research, we introduced a Bayesian statistical test based on Kibble's bivariate gamma distribution (BKD) [9] to tackle the detection challenge under consideration. This method responded to the need to adjust confidence in the stationarity of the radiological background signal by integrating absolute and relative features from two Bayesian tests [10] developed for the same application. However, the representativeness of the BKD test is constrained by its reliance on comparing two discrete distributions, assuming independence of counting spectra channel contents. In contrast, while the classical frequentist Kolmogorov-Smirnov (KS) test performs less optimally for our detection problem [11], it does not assume independence or interchangeability of histogram bins but operates with samples from continuous underlying distributions. In this study, our aim is to formulate a statistical test that considers the correlation between samples arranged in a histogram. This test aims to offer higher specificity and sensitivity than both the KS test and its discrete counterpart (BKD), addressing the requirements and limitations of large-scale decommissioning projects. Thus, we developed a new Bayesian hypothesis test based on a Pólya tree-type *a priori*, referred to as the Bayesian test by mixture of Kibble-Pólya Tree (BKPT). The innovation lies in the adoption

of a finite Pólya tree model, which introduces correlations between the relative contents of spectral channels based on their positions. This approach is particularly suitable for high-resolution spectrometers like HPGe detectors, assuming full energy peaks (FEP) exhibit a shape closely resembling a Dirac distribution. However, this assumption becomes less robust for lower-resolution detectors like Sodium iodide (NaI(Tl)) or Lanthanum bromide (LaBr₃) scintillators, which offer greater scalability and field adaptability. The BKPT test, being a more comprehensive Bayesian test, enables the processing of spectra where the correlation between adjacent counting channels' contents cannot be ignored. This is particularly relevant for lower-resolution detectors, as detailed in Section 3. It accommodates output spectra from both high- and low-resolution detections, under the assumption that the spectrometer used in this study is the same as the HPGe detector described in the features of the BKD test [10]. Implementing the BKD and BKPT tests involves modeling not only the null hypothesis, labeled, $H_0$, for no uranium contamination, as in frequentist tests, but also the alternative hypothesis, labeled, $H_1$, based on foreknowledge about the expected contamination signal. This prediction relies on prior acquisitions within the nuclear facility to be decommissioned and MCNP6 simulation outputs [12]. Establishing precise criteria for rejecting the null hypothesis and triggering a uranium contamination alarm is crucial in Bayesian decision-making, relying on the posterior probabilities of the null and alternative hypotheses. We evaluate the specificity and sensitivity of the BKPT test by considering a gamma-ray source term model for various uranium surface contaminations, and we compare its performance with the BKD test using standard receiver operating characteristics (ROC) curves [13].

There are five main sections in this paper: in Section 2, we briefly present the physical model of our detection problem; in Section 3, we introduce the mathematical formulations of the BKD and BKPT tests; in Section 4, we discuss the respective specificities and sensitivities of both tests, based on findings from a numerical study fed with experimental data.

## Physical model for detecting Low-Activity uranium contamination

In this section, we provide a brief overview of the key components of the physical model incorporated into the BKPT test.

### Uranium gamma-ray peaks contamination

As previously stated, our approach to detecting low-activity uranium contamination involves using high-resolution gamma-ray spectrometry, specifically targeting the most intense emission lines from the decay chains of U-235 and U-238 [14]. The selection of these lines is based on the detector's full energy peak efficiency at the associated energies and the presence of neighboring FEP from the radiological background. To capture the counting signal, we utilized an n-type coaxial HPGe detector [15]. This HPGe detector is equipped with a 2-mm aluminum window to safeguard the crystal from potential surface irregularities.

Under these specified conditions, and relying on previous on-site acquisitions [9,10,11], we narrowed down the list of gamma-ray peaks of interest to a primary peak at 185.7 keV and three additional peaks at {143.8: 163.4; 205.3} keV for the U-235 chain, along with a main emitter peak at 1001.0 keV for the U-238 chain. Fig. 1 illustrates a graphical representation of the U-235 and U-238 peaks of interest.

As extensively discussed in reference [10], the process of selecting Regions of Interest (ROI) within the recorded spectra, ranging from 0 to 3 MeV, was conducted following the energy and resolution calibration of the spectrometry chain. This calibration was performed using a reference europium-152 (Eu-152) laboratory source. The finalized list of spectral ROIs pertinent to the uranium gamma-ray signature is as follows: $ROI_1 = [140-148]$ keV, $ROI_2 = [160-166]$ keV, $ROI_3 = [182-188]$ keV, $ROI_4 = [202-208]$ keV, and $ROI_5 = [998-1004]$ keV.
<figure><img src="image_3.png" /><figcaption>Fig. 1. Cumulative spectra linked to the background noise and uranium samples with enrichment levels of 0.7% and 9.5%. (a) [50–250] keV, (b) [750–1100] keV.</figcaption></figure>

The BKD test exhibits the capability to handle discontinuous spectral ROIs, a characteristic that will be clarified in Section 3. However, the underlying rationale behind the BKPT test, which involves accounting for counting correlations between neighboring channels, requires operating over a continuous region. Therefore, we opted for the continuous range $ROI = [140–208]$ keV, which includes the four prominent lines from U-235 decay, to evaluate the performance of the newly introduced BKPT test. To ensure a fair comparison, we constrained the spectral region for assessing the BKD to the discontinuous union $ROI_1 \cup ROI_2 \cup ROI_3 \cup ROI_4$.

### Detector response estimation using MCNP simulation

The three following steps outline the procedure we applied to compute the expected response of the HPGe detector to the emission lines of interest. Firstly, we built a transfer function relating the uranium minimum surface activities to be detected (MSAD, in Bq.m$^{-2}$ over $2\pi$) to the emission rates (in s$^{-1}$) of the gamma-ray lines identified in Section 2.1 [11]. This function was constructed assuming the existence of secular equilibrium between U-238 and Pa-234 m, which generates a 1001.0 keV gamma-ray upon decay. The transfer function depends on two parameters: the intensity of the gamma-ray line and the enrichment

level of U-235, expressed as mass fraction (wt%). In alignment with reference [11], the MSAD was set equal to the most restrictive decommissioning criterion applicable to the site where the counting data were acquired, specifically, 500 Bq.0.1 m$^{-2}$ over $2\pi$, as defined by the French Nuclear Safety Authority. This ensures the detection sensitivity meets regulatory requirements. Using the parameters of gamma-ray line intensity and U-235 enrichment level, we applied the transfer function to calculate the expected response of the HPGe detector. This step ensures that the detector's response is accurately modeled based on realistic environmental and operational conditions.

Secondly, we constructed a numerical model of the HPGe detector using MCNP6.1 [12] and calibrated it with uranium test samples of known activities and U-235 enrichment levels. This calibration allowed us to estimate count rates (in s$^{-1}$) within predefined ROIs for a spread contamination of 500 Bq.0.1 m$^{-2}$. Fig. 2 represents the orientation and structure of the simulated HPGe detector. The HPGe detector is shown in the lower right corner, while the model photon source is evenly distributed across the 0.1 m$^2$ concrete block surface in the upper left corner. The photon flux, derived from the MSAD using the enrichment-dependent transfer function, generates gamma-ray fluxes uniformly over a 0.1 m$^2$ surface at a distance of $d = 22.36$ cm from the detector's

<figure><img src="image_4.png" /><figcaption>Fig. 2. Simulated HPGe detector structure using MCNP6.1 (a) and its orientation (b).</figcaption></figure>
active area. Under these conditions, the response to the main 185.7 keV line is expected to be on the order of $10^{-2}$ s$^{-1}$.

The final step involves selecting relevant U-235 enrichments corresponding to three subunits (SU) of the nuclear installation to be decommissioned: SU1, SU2, and SU3. Average expected U-235-enrichment levels of 0.62 wt%, 2.03 wt% and 6.84 wt% were associated with these subunits, respectively. These expected enrichments will serve as the basis for generating a priori spectral signatures from uranium contaminations for the implementation of Bayesian tests. Additionally, the FEP responses of the detector to U-235 mass fractions of 0.7 wt%, 3 wt% and 8 wt% were estimated to generate tested spectra for simulating real-time, random measurements in SU1, SU2, and SU3, as detailed in Section 2.4. It should be noted that the signatures of Radium-226, Bismuth-214, and Plomb-214 (Ra-226, Bi-214, and Pb-214) were not considered in this study because they are associated with descendants of U-238, whereas our primary focus remains on detecting U-235. Therefore, we investigated various levels of potential U-235 enrichment. The chosen lines are sufficient for achieving satisfactory results, but we may add more in the future. Additionally, we used the same simulated signal as in previous studies to ensure a fair comparison with the BKD test. Fig. 3 presents a diagram outlining the three steps of estimating detector response through MCNP simulation.

## Background radiation variation

Detecting a uranium signal within an enclosed concrete environment is particularly challenging due to the presence of a radiological background. This background primarily arises from potassium 40 (K-40) and the decay chains of natural uranium-238 (U-238), uranium-235 (U-235), and thorium-232 (Th-232) [16]. These signatures characterize mainly the gamma-ray background, as determined by in-situ measurements conducted within the nuclear facility. The background appears as broad diffusion continua or adjacent FEP, complicating the detection process. The challenge is further amplified by variations in the intensity or normalized distribution of this background across different facility locations. As detailed in reference [10], our initial task was to quantify this variability across the three subunits SU to be decommissioned.

To achieve this, we computed relative deviations between spectra recorded at different spots within each SU, maintaining a consistent detector distance ($d = 22.36$ cm) from every investigated $0.1\text{m}^2$ square surface. This computation allowed us to characterize variations in both the amplitude and shape of the radiological background, considering

non-normalized and normalized counting spectra, respectively. Our analysis of the collected data indicated that, on average, the background amplitude and shape over the specified ROIs are expected to vary by approximately 1% across one SU and locally up to 10%. Consequently, the procedure for assessing the performance of Bayesian tests in this application will need to model fluctuations within these orders of magnitude.

## Random generation of tested spectra

Our objective in implementing Bayesian detection tests is to determine, based on an integrated count spectrum measurement, whether the detector response is solely due to the radiological background or if it also includes an isotopic signal indicative of a uranium surface contamination. To evaluate the effectiveness of these tests, we generated random test spectra that simulate the spectrometer's response in both scenarios [10].

The numerical process employed time-independent probability histograms, allowing the acquisition time $t$ to be treated as a variable parameter in the simulator. The process for each subunit SU to be decommissioned is outlined as follows [10]:

j) One of the background spectra acquired inside the SU is normalized to serve as uranium-free, reference spectrum;

ii) The expected spectral signature of a uranium surface activity equal to MSAD, with the average U-235 enrichment known for the SU (0.62 wt%, 2.03 wt%, and 6.84 wt% for SU1, SU2, and SU3 respectively), serves as a net, a priori contamination spectrum;

iii) The time-independent probability histogram associated with a test spectrum with no uranium contamination signal is formed by the random and uniformly sampled combination of two background spectra acquired inside the SU;

iv) The time-independent probability histogram associated with a test spectrum containing an additional uranium contamination signal is formed by the sum of a distribution ii) and the expected spectral signature of a uranium surface activity equal to MSAD with a U-235 enrichment that differs from the average enrichment level known for SU1, SU2, and SU3 respectively;

v) The construction of uranium-free spectra (to estimate the false alarm rate (FAR) of the procedure) and uranium-rich spectra (to estimate the true detection rate (TDR) of the procedure) is performed by multiplying one of histograms iii) or iv) by a preset integration time $t$, followed by a random draw according to a Poisson distribution within the thus-constructed, time-dependent expectation vector.

## Statistical detection tests based on Bayesian inference

Like in any hypothesis testing, a detection procedure using Bayesian inference aims to decide, based on the counting spectrum recorded over an integration time labeled $t$, whether:

- the observer is in the presence of solely the radiological background (a scenario corresponding to the null hypothesis, denoted as $H_0$);

- or in the presence of both the radiological background and an additional isotopic signal to be detected (referred to as the alternative hypothesis, $H_1$).

In this section, we present the theory and mathematical formulation of the newly introduced BKPT test, following a review of the essential features of the BKD test [9], which will serve as a comparative detection procedure.

### Bayesian test by mixture of Kibble-Dirichlet (BKD)

The construction of the BKD test requires some foreknowledge of

<figure><img src="image_5.png" /><figcaption>Fig. 3. Diagram of the detector response estimation steps using MCNP simulation.</figcaption></figure>
both spectral signatures of the radiological background and the expected uranium contamination inside a given SU. The former is obtained through on-site acquisitions carried out over a duration per spot labelled $T_{ref}$, while the latter is derived from Monte Carlo simulations [10]. The keystone of the BKD test lies in the statistical level introduced between these a priori spectral signatures and the measurements. It treats the count histograms acquired over the reference acquisition time $T_{ref}$ and the test acquisition time $T_{test}$ as outcomes of sampling from multivariate Poisson laws with random intensities. These intensities are generated according to the product of a Dirichlet variable and a positive scale factor, eventually modeled by a bivariate gamma distribution of Kibble [9]. The notations and statistical models underlying the formulation of the BKD test are as follows:

(i) Multichannel countings $m_{ref}$ and $m_{test}$ are respectively associated with the acquisitions carried out over $T_{ref}$ and $T_{test}$. We note $m_{ref} = (m_{ref}[1], \dots, m_{ref}[K])$ and $m_{test} = (m_{test}[1], \dots, m_{test}[K])$, where $m_{ref}[k]$ and $m_{test}[k]$ respectively denote the count numbers recorded in the $k^{th}$ channels of the reference and tested spectra, and K represents the number of exploited channels in these spectra.

(ii) The number of counts in the $k^{th}$ channel of a histogram originating only from the radiological background is modelled as random draws from Poisson distributions with respective vector parameters $\lambda_{Bref} \cdot x_B[k]$ and $\lambda_{Btest} \cdot x_B[k]$ for the reference and test spectra, respectively. The vector $x_B$ is supposed to take values in a Dirichlet distribution whose vector parameter is given by the product $\gamma_B = K \cdot \nu_B \cdot s_B$, where vector $s_B$ denotes the normalized a priori spectral signature of the background, and the scalar $\nu_B$ is a weighting factor reflecting the user's confidence in the a priori information.

(iii) A key aspect of the BKD test is that the counting intensities $\lambda_{Bref}$ and $\lambda_{Btest}$ are postulated to follow a Kibble bivariate gamma distribution, referred to as Kibb, according to:

$$ \lambda_{Bref}, \lambda_{Btest} \sim \text{Kibb}(a_B, b_B, \rho) \quad (1) $$

and the probability density function, denoted $f_{Kibb}$, is explained [17,18] as:

$$ f_{\text{Kibb}}(\lambda_{\text{Bref}}, \lambda_{\text{Btest}}) = \frac{b_B^{2a_B}}{\Gamma(a_B)(1-\rho)^{a_B}} \times \exp\left(-\frac{b_B(\lambda_{\text{Bref}} + \lambda_{\text{Btest}})}{1-\rho}\right) \times \sum_{i=0}^{\infty} \frac{(\lambda_{\text{Bref}} \lambda_{\text{Btest}})^{a_B+i-1}}{i! \Gamma(a_B+i)} \times \rho^i \left(\frac{b_B}{1-\rho}\right)^{2i} \quad (2) $$

correlation coefficient $\rho$ is defined by the user and expresses the a priori correlation between radiological background intensities over the reference period $T_{ref}$ and the measurement period $T_{test}$ (i.e., the stationarity of the background). Note that $\rho = 0$ corresponds to assuming that such intensities are uncorrelated, meaning $\lambda_{Bref}$ and $\lambda_{Btest}$ are two independent gamma parameters. Conversely, taking $\rho \to 1$ provides a stationarity postulate: $\lambda_{Bref} \approx \lambda_{Btest}$ [9].

(v) In a similar fashion, the number of counts in the $k^{th}$ channel of a histogram generated solely by a uranium contamination is assumed to follow a Poisson distribution with vector parameter $\lambda_S \cdot x_S[k]$. The counting intensity $\lambda_S$ from a uranium signal in the measured spectrum $m_{test}$ is postulated to take values in a gamma distribution with parameters $(a_S; b_S)$, where $a_S$ is set equal to $1/2$ and $b_S$ is inversely proportional to an a priori average intensity of the contamination signal. The vector $x_S$, on the other hand, follows a Dirichlet distribution whose vector parameter is given by the product $\gamma_S = K \cdot \nu_S \cdot s_S$. The vector $s_S$ denotes the normalized a priori spectral signature of the uranium signal, and the scalar $\nu_S$ is a weighting factor reflecting the degree of confidence of the user in said a priori.

(vi) The Bayesian inference procedure aims to determine the a posteriori probabilities that the investigated surface has the same radiological status as the reference surface or to presents a uranium contamination, recalling that:

- Under hypothesis $H_0$, both spectra $m_{ref}$ and $m_{test}$ consist solely of counts from the radiological background;

- Under hypothesis $H_1$, the spectrum $m_{ref}$ contains only radiological background counts, whereas the spectrum $m_{test}$ includes contributions from both the radiological background and a uranium contamination signal.

Calculating the joint conditional probabilities $\mathbb{P}(m_{ref}, m_{test}|H_0)$ and $\mathbb{P}(m_{ref}, m_{test}|H_1)$ is a prerequisite for accepting either $H_0$ or $H_1$. To achieve this, we compute the joint likelihood of $m_{ref}$, $m_{test}$, $x_B$, $\lambda_{Bref}$, and $\lambda_{Btest}$ under $H_0$, as well as the joint likelihood of $m_{ref}$, $m_{test}$, $x_B$, $\lambda_{Bref}$, $\lambda_{Btest}$, $x_S$, and $\lambda_S$ under $H_1$.

The marginalization of the first joint conditional likelihood directly yields the following expression for $\mathbb{P}(m_{ref}, m_{test}|H_0)$, where $\mathbb{F}_1$ denotes the ordinary hypergeometric function [9]:

$$ \mathbb{P}(\mathbf{m}_{ref}, \mathbf{m}_{test}|H_0) = \frac{T_{ref}^{M_{ref}} T_{test}^{M_{test}}}{\prod_{k=1}^{K} m_{ref}[k]! m_{test}[k]! \Gamma(\mathbf{M}_{ref} + \mathbf{M}_{test} + \mathbf{K}_{VB})} \frac{\Gamma(\mathbf{K}_{VB})}{\prod_{k=1}^{K} \frac{\Gamma(m_{ref}[k] + m_{test}[k] + \gamma_B[K])}{\Gamma(\gamma_B[k])}} \times \left(\frac{b_B^{a_B}}{\Gamma(a_B)}\right)^2 \frac{\Gamma(a_B + \mathbf{M}_{ref}) \Gamma(a_B + \mathbf{M}_{test})}{(1-\rho)^{a_B} \tilde{T}_{ref}^{a_B + \mathbf{M}_{ref}} \tilde{T}_{test}^{a_B + \mathbf{M}_{test}}} \times (1-z)^{-(a_B + \mathbf{M}_{ref} + \mathbf{M}_{test})} \mathbb{F}_1(-\mathbf{M}_{ref}, -\mathbf{M}_{test}, a_B, z) \quad (4) $$

where $\Gamma$ refers to the gamma function.

(iv) From the Eqs. (1) and (2), the first scalar variable $a_B > 0$ is defined as $1/2$, whereas the second variable $b_B > 0$ is inversely proportional to the a priori average intensity of the background. The last variable $\rho \in [0; 1]$ denotes the Pearson correlation coefficient, representing the correlation between $\lambda_{Bref}$ and $\lambda_{Btest}$, and can be expressed as

$$ \rho = \frac{\mathbb{E}(\lambda_{Bref} \lambda_{Btest}) - \bar{\lambda}_{Bref} \bar{\lambda}_{Btest}}{\sigma_{Bref} \sigma_{Btest}} \quad (3) $$

where $\bar{\lambda}_{Bref}$ and $\sigma_{Bref}$ represent the mean and standard deviation of the variable $\lambda_{Bref}$, $\bar{\lambda}_{Btest}$ and $\sigma_{Btest}$ represent the mean and standard deviation of the variable $\lambda_{Btest}$, and $\mathbb{E}$ denotes the expectation operator. The

where $M_{ref} = \sum_{k=1}^{K} m_{ref}[k]$, $M_{test} = \sum_{k=1}^{K} m_{test}[k]$, $\tilde{T}_{ref} = T_{ref} + \frac{b_B}{1-\rho}$, $\tilde{T}_{test} = T_{test} + \frac{b_B}{1-\rho}$, and $z = \frac{\rho}{T_{ref}T_{test}} \left(\frac{b_B}{1-\rho}\right)^2$.

The marginalization of the second joint conditional likelihood with both $x_B$ and $x_S$ is not as straightforward as the previous one and requires the integration and development of a multivariate polynomial. This polynomial takes advantage of the moment structure of the Dirichlet distribution, which allows a coefficient-by-coefficient identification of the powers of the polynomial to various instantiations of the gamma function [9]. Eventually, an expression for $\mathbb{P}(m_{ref}, m_{test}|H_1)$ may be given by:
$$
\begin{equation}
\begin{aligned}
\mathbb{P}(\mathbf{m}_{\text{ref}}, \mathbf{m}_{\text{test}}|H_1) = & \frac{T_{\text{ref}}^{M_{\text{ref}}} T_{\text{test}}^{M_{\text{test}}}}{\prod_{k=1}^{K} m_{\text{ref}}[k]! m_{\text{test}}[k]! \prod_{k=1}^{K} \Gamma(\gamma_B[k]) \prod_{k=1}^{K} \Gamma(\gamma_S[k])} \frac{1}{\Gamma(a_S)} \left(\frac{b_S}{T_{\text{test}} + b_S}\right)^{a_S} \left(\frac{b_B^{a_B}}{\Gamma(a_B)}\right)^2 \\
& \times \frac{\Gamma(a_B + \mathbf{M}_{\text{ref}})}{(1 - \rho)^{a_B} \tilde{T}_{\text{ref}}^{a_B + M_{\text{ref}}} \tilde{T}_{\text{test}}^{a_B + M_{\text{test}}}} (1 - z)^{-(a_B + M_{\text{ref}} + M_{\text{test}})} \\
& \times \sum_{j=0}^{M_{\text{test}}} \frac{\omega_j \Gamma(j + a_S) \Gamma(M_{\text{test}} - j + a_B)_2 F_1(-M_{\text{ref}}, -M_{\text{test}} + j, a_B, z) \psi^j}{\Gamma(\kappa_{\text{V}_S} + j) \Gamma(\kappa_{\text{V}_B} + \mathbf{M}_{\text{ref}} + \mathbf{M}_{\text{test}} - j)}
\end{aligned}
\tag{5}
\end{equation}
$$

where $M_{\text{ref}}$, $M_{\text{test}}$, $\tilde{T}_{\text{ref}}$, $\tilde{T}_{\text{test}}$, and $z$ are defined as before, and $\psi = \frac{(1-z)\tilde{T}_{\text{test}}}{\tilde{T}_{\text{test}}+b_S}$.

The expression of $\mathbb{P}(\mathbf{m}_{\text{ref}}, \mathbf{m}_{\text{test}}|H_1)$ uses coefficients $\omega_j$ that are derived in two steps. First, we form a set of $K$ vectors labeled $(\mathbf{V}_k)_{1 \le k \le K}$, where the $j^{\text{th}}$ coordinate $V_k[j]$ of each of these vectors is calculated as:

$$ V_k[j] = \begin{pmatrix} m_{\text{test}}[k] \\ j \end{pmatrix} \frac{\Gamma(j + \gamma_S[k])}{\Gamma(\gamma_S[k])} \frac{\Gamma(m_{\text{ref}}[k] + m_{\text{test}}[k] - j + \gamma_B[k])}{\Gamma(m_{\text{ref}}[k] + m_{\text{test}}[k] + \gamma_B[k])} \quad (6) $$

Second, the vector $\Omega = (\omega_j)_{1 \le j \le M_{\text{test}}}$ is obtained by $K$ discrete convolutions of the vectors $(\mathbf{V}_k)_{1 \le k \le K}$, so that: $\Omega = \mathbf{V}_1 \otimes \cdots \otimes \mathbf{V}_k \cdots \otimes \mathbf{V}_K$. The convolution product symbolized in the preceding expression is made explicit as:

$$ \forall (m, n) \in [1, K]^2, \forall q \ge 1, (V_m \otimes V_n)[q] = \sum_{q=-\infty}^{+\infty} V_m[q] \bullet V_n[p-q] \quad (7) $$

The ratio of the a posteriori probabilities for both the null and alternative hypotheses is then obtained by applying Bayes' theorem:

$$ \frac{\mathbb{P}(H_0 | \mathbf{m}_{\text{ref}}, \mathbf{m}_{\text{test}})}{\mathbb{P}(H_1 | \mathbf{m}_{\text{ref}}, \mathbf{m}_{\text{test}})} = \frac{\mathbb{P}(\mathbf{m}_{\text{ref}}, \mathbf{m}_{\text{test}} | H_0) \mathbb{P}(H_0)}{\mathbb{P}(\mathbf{m}_{\text{ref}}, \mathbf{m}_{\text{test}} | H_1) \mathbb{P}(H_1)} \quad (8) $$

The so-called prior probabilities $\mathbb{P}(H_0)$ and $\mathbb{P}(H_1)$ are considered as equal in this study. Additionally, we have $\mathbb{P}(H_1|\mathbf{m}_{\text{ref}}, \mathbf{m}_{\text{test}}) = 1 - \mathbb{P}(H_0|\mathbf{m}_{\text{ref}}, \mathbf{m}_{\text{test}})$.

From these equalities and Eq. (8), we obtain the following key expression for the a posteriori probability of the null hypothesis:

$$ \mathbb{P}(H_0 | \mathbf{m}_{\text{ref}}, \mathbf{m}_{\text{test}}) = \frac{1}{1 + \kappa_{\text{BKD}}} \quad (9) $$

where $\kappa_{\text{BKD}}$ is the Bayes factor associated with the BKD test, formally defined as:

$$ \kappa_{\text{BKD}} = \frac{\mathbb{P}(\mathbf{m}_{\text{ref}}, \mathbf{m}_{\text{test}}|H_1)}{\mathbb{P}(\mathbf{m}_{\text{ref}}, \mathbf{m}_{\text{test}}|H_0)} \quad (10) $$

and explicitly derived from Eqs. (4) and (5) as:

$$ \kappa_{\text{BKD}} = \frac{\Gamma(\kappa_{\text{V}_S})}{\prod_{k=1}^{K} \Gamma(\gamma_S[k])} \frac{\Gamma(\mathbf{M}_{\text{ref}} + \mathbf{M}_{\text{test}} + \kappa_{\text{V}_B})}{\prod_{k=1}^{K} (m_{\text{ref}}[k] + m_{\text{test}}[k] + \gamma_B[k])} \left(\frac{b_S}{T_{\text{test}} + b_S}\right)^{a_S} \\
\times \frac{1}{\Gamma(a_S)\Gamma(a_B + \mathbf{M}_{\text{test}})_2 F_1(-\mathbf{M}_{\text{ref}}, -\mathbf{M}_{\text{test}}, a_B, z)} \\
\times \sum_{j=0}^{M_{\text{test}}} \frac{\omega_j \Gamma(j + a_S) \Gamma(\mathbf{M}_{\text{test}} - j + a_B)_2 F_1(-\mathbf{M}_{\text{ref}}, -\mathbf{M}_{\text{test}} + j, a_B, z) \psi^j}{\Gamma(\kappa_{\text{V}_S} + j) \Gamma(\kappa_{\text{V}_B} + \mathbf{M}_{\text{ref}} + \mathbf{M}_{\text{test}} - j)} \quad (11) $$

Eventually, the BKD hypothesis test against a given $\alpha$ - risk (standing for the false alarm rate, abbreviated FAR), reads as:

- If $\frac{1}{1+\kappa_{\text{BKD}}} < \eta_a$, then $H_0$ is rejected and a uranium contamination signal is detected;

- Otherwise, $H_0$ is accepted.

We note that the relation between the decision threshold $\eta_a$ and the

$\alpha$ - risk is determined numerically via a Monte Carlo method.

## Bayesian test by mixture of Kibble-Pólya tree (BKPT)

In the BKD test, the assumption of independence among the contents of spectral channels in multichannel countings $\mathbf{m}_{\text{ref}}$ and $\mathbf{m}_{\text{test}}$ is based on the underlying probability distributions. These counts follow a Poisson distribution, with parameters governed by a discrete Dirichlet distribution. This approach is suitable for high-resolution spectrometers like HPGe detectors, where each FEP approximates a Dirac distribution. However, this assumption becomes problematic when dealing with lower-resolution detectors such as (NaI(Tl)) or LaBr<sub>3</sub> scintillators, which may offer greater scalability and field adaptability but may not exhibit FEP shapes resembling Dirac distributions.

The motivation behind introducing a more universal Bayesian test is to accommodate spectra where the correlation between neighboring counting channels cannot be ignored. Accounting for this interdependence may lead to improved figures of merit (FOM), particularly for high-resolution spectrometer such as the HPGe detector presented in Section 2. As we'll explore in this section, using a Pólya tree-type model instead of the Dirichlet model in the BKD test can address these expectations.

Pólya trees form a class of distributions commonly used as a priors in nonparametric Bayesian modeling, particularly in scenarios where random probabilities need to be measured [19,20]. A finite Pólya tree (FPT), shares similarities with the Dirichlet distribution in that both relate to a $k$ - simplex, denotes as $\mathbf{x} = (x_k)_{1 \le k \le K}$, representing either a normalized radiological background or a uranium signal spectrum. However, the crucial distinction lies in the FPT's ability to establish correlations among the relative contents of spectral channels based on their positions.

One approach to achieve this is by assuming that the support of the counting spectrum $\mathbf{x}$ is obtained from repeated and random dyadic cuts. Consequently, constructing a binary FPT requires that the number of channels $K$ forming a continuous ROI be a power of 2: $K = 2^L$, with $L$ denotes the number of levels (i.e., cutting sequences) of the FPT. For each level $l$ of the tree, the number of calculated proportions, labeled $K_l$, will also be a power of 2: $\forall l \in [1, L]$, $K_l = 2^l$, with $K_L = K$ by definition of $L$.

For the sake of simplicity, we rename the final proportions in any a priori normalized spectrum $\mathbf{x}$ from 0 to $K-1$. The model assumes that a normalized counting value $x_k$ is given by the product of all relative proportions $p_{l,\Delta(L-l,k)}$, where $1 \le l \le L$ and $\Delta(L-l, k) = |2^{-(L-l)}K|$, representing the initial number of counts $p_{0,0} = 1$ retained through the sequence of cuts leading to the formation of this channel:

$$ \forall k \in [0, K-1], x_k = \prod_{l=1}^{L} p_{l, \Delta(L-l, k)} \quad (12) $$

<caption>Table 1</caption>

3 - level FPT modeling an 8 - channel normalized spectrum.

<table><tr><td colspan="4">p<sub>0,0</sub> = p<sub>1,0</sub> + p<sub>1,1</sub> = 1</td></tr><tr><td colspan="2">p<sub>1,0</sub></td><td colspan="2">p<sub>1,1</sub></td></tr><tr><td>p<sub>2,0</sub></td><td>p<sub>2,1</sub></td><td>p<sub>2,2</sub></td><td>p<sub>2,3</sub></td></tr><tr><td>p<sub>3,0</sub></td><td>p<sub>3,1</sub></td><td>p<sub>3,2</sub></td><td>p<sub>3,3</sub></td><td>p<sub>3,4</sub></td><td>p<sub>3,5</sub></td><td>p<sub>3,6</sub></td><td>p<sub>3,7</sub></td></tr><tr><td>x<sub>0</sub></td><td>x<sub>1</sub></td><td>x<sub>2</sub></td><td>x<sub>3</sub></td><td>x<sub>4</sub></td><td>x<sub>5</sub></td><td>x<sub>6</sub></td><td>x<sub>7</sub></td></tr></table>
**Table 1** illustrates the modeling of an 8 – Channel normalized spectral ROI **x** = (x<sub>k</sub>)<sub>0≤k≤7</sub> by a 3 – level FPT.

The keystone of the model is the statistical distributions governing the sequence of dyadic cuts, represented by the 2K − 2 strictly positive coefficients (p<sub>l,2k</sub>, p<sub>l,2k+1</sub>)<sub>1≤l≤L,0≤k≤K−1</sub>. The FPT model assumes that the content of any bin the normalized energy histogram is obtained by the product of samples from L independent random variable following a beta distribution, labeled B:

$$
\forall l \in [1, L], \forall k \in [0, K-1], \quad \begin{array}{l} p_{l,2k} \sim \text{B}(\cdot|s_{l,2k}, s_{l,2k+1}) \\ p_{l,2k+1} = 1 - p_{l,2k} \end{array} \qquad (13)
$$

In Eq. (13), the recursive construction of parameters such as s<sub>l,2k</sub> and s<sub>l,2k+1</sub> will be explained further in the following subsection. It's notable from this equation that a Dirichlet distribution is essentially a FPT [19,20], where for every *l* and *k*, the equality holds: s<sub>l,k</sub> = s<sub>l+1,2k</sub> + s<sub>l+1,2k+1</sub>. This reduction in the number of free parameters to *K* is why the FPT is considered a generalization of the Dirichlet distribution. It's worth emphasizing that Pólya trees, as prior distributions on density functions, can be viewed as FPTs with an infinite number of levels L→∞.

Construction of the BKPT test

In the context of the BKPT test, the null and alternative hypotheses H₀ and H₁ maintain the same definitions as outlined in Section 3.1. All relevant parameters and data retain the same notations as before. The statistical assumptions established during the construction of the BKD test remain intact in the formulation of the BKPT test, with notable exception: the distributions governing the vectors **x**<sub>B</sub> and **x**<sub>S</sub>. Specifically, while the counting intensities λ<sub>B,ref</sub> and λ<sub>B,est</sub> continue to be assumed to follow a bivariate gamma distribution of Kibble, the key difference lies in the distributional assumptions for the vectors **x**<sub>B</sub> and **x**<sub>S</sub>. Unlike the BKD test, in the BKPT test, these vectors are no longer presumed to be distributed according to Dirichlet distributions but rather finite Pólya trees, labeled FPTs.

As previously mentioned, a FPT is a probability distribution defined over a strictly positive, normalized simplex **x**. Like any distribution, it requires parameters for its construction. Let's denote **S**<sub>B</sub> as the vector of parameters governing the normalized distribution of the radiological background, FPT(**x**<sub>B</sub> | **S**<sub>B</sub>). The only constraint when constructing the parameters **S**<sub>B</sub> = (b<sub>l,k</sub>)<sub>1≤l≤L,0≤k≤K−1</sub>, with K = 2<sup>L</sup>, is that the mean of FPT(**x**<sub>B</sub> | **S**<sub>B</sub>) equals the *a priori* spectral signature of the background, labeled as **s**<sub>B</sub> in Section 3.1.

To facilitate the construction of the vector S<sub>B</sub>, we introduce two strictly positive scale parameters σ and c<sub>0</sub>. A higher value of σ results in a more rigid (in terms of low-variance) FPT. As previously mentioned, FPT(x<sub>B</sub> | S<sub>B</sub>) needs to be parametrized by 2K − 2 strictly positive coefficients. Considering K = 2<sup>L</sup>, assigning individual values to every (b<sub>l,k</sub>)<sub>0≤l≤L,0≤k≤K−1</sub> coefficient is impractical. However, it can be shown that as long as:

$$
\forall l \in [1, L], \forall k \in [0, K-1], b_{l-1,k} \propto b_{l,2k} + b_{l,2k+1} \quad (14)
$$

where, for every level of the tree *l* ∈ [1, *L*], *K*<sub>*l*−1</sub> = 2<sup>*l*−1</sup> labels the num-
ber of calculated proportions. Then, the expected value or mean *E*(**x**<sub>B</sub>) of
the FPT corresponds to the normalized coefficients of the *L*<sup>th</sup> level, so:

$$
\mathbb{E}(\mathbf{x}_B) = \frac{b_{l...}}{C} = \frac{\{b_{L,k}\}_{0 \le k \le K-1}}{C} \qquad (15)
$$

with $C = \sum_{k=0}^{K-1} b_{L,k}$.

The construction rule that this mean be equal to the *a priori* signature
of the radiological background then translates as:

$$
\forall k \in [0, K-1], \frac{b_{L,k}}{C} = s_B[k] \Leftrightarrow b_{L,k} = C \cdot s_B[k] \quad (16)
$$

The simplest way to abide by Eq. (14) is to introduce a unique

multiplication factor to move from any level *l*<sup>th</sup> down to the (l − 1)<sup>th</sup> level of the tree. The incremental scale factor, which we labeled σ, will play this part for every level of the FPT, so that:

$$
\forall l \in [1, L], \forall k \in [0, K_{l-1}], b_{l-1,k} = \frac{b_{l,2k} + b_{l,2k+1}}{2\sigma} \quad (17)
$$

Eventually, we set the first dyadic cutting coefficient to the initial scale parameter, which we labeled c₀. It follows from Eqs. (16) and (17) that the normalization constant.

$$
C = c_0 \cdot (2\sigma)^L \qquad (18)
$$

In a nutshell, the FPT parameters are then constructed as follows:

b<sub>0,0</sub> = c<sub>0</sub>; ∀k ∈ [0, K − 1], b<sub>L,k</sub> = s<sub>B</sub>[k] • c<sub>0</sub> • (2σ)<sup>L</sup>. Then, recursively, from bottom to top of the tree:

$$
\forall l \in [1, L], \forall k \in [0, K_{l-1}], b_{l-1,k} = \frac{b_{l,2k} + b_{l,2k+1}}{2\sigma} \quad (19)
$$

The finite Pólya tree modeling the normalized distribution of a uranium contamination signal FPT(**x**<sub>s</sub> | **S**<sub>s</sub>) is constructed in the same manner. Once we require that the mean of the FPT be equal to the *a priori* signature of a uranium contamination signal within the investigated SU, labeled **s**<sub>s</sub>, the vector of parameters **S**<sub>s</sub> = (s<sub>l,k</sub>)<sub>1≤l≤L,0≤k≤K−1</sub>, with K = 2<sup>L</sup>, is thus built as follows:

s<sub>0,0</sub> = c<sub>0</sub>; ∀k ∈ [0, K − 1], s<sub>L,k</sub> = s<sub>S</sub>[k] • c<sub>0</sub> • (2σ)<sup>L</sup>. Then, recursively, from bottom to top of the tree:

$$
\forall l \in [1, L], \forall k \in [0, K_{l-1}], s_{l-1,k} = \frac{s_{l,2k} + s_{l,2k+1}}{2\sigma} \quad (20)
$$

The formulation of the Bayesian hypothesis test, similar to the one for BKD in Section 3.1, requires the calculation of both joint conditional probabilities P(m<sub>ref</sub>, m<sub>test</sub> | H₀) and P(m<sub>ref</sub>, m<sub>test</sub> | H₁).

First, the joint likelihood of **m**<sub>ref</sub>, **m**<sub>test</sub>, **x**<sub>B</sub>, **λ**<sub>B,ref</sub>, and **λ**<sub>B,est</sub> under the
null hypothesis H<sub>0</sub> is expressed as:

$$
\begin{equation}
\begin{split}
f(\mathbf{m}_{\text{ref}}, \mathbf{m}_{\text{test}}, \lambda_{\text{B,ref}}, \lambda_{\text{B,est}}, \mathbf{x}_{\text{B}} | H_0) ={}& \frac{T_{\text{ref}}^{M_{\text{ref}}} T_{\text{test}}^{M_{\text{test}}}}{\prod_{k=1}^{K} m_{\text{ref}}[k]! m_{\text{test}}[k]!} e^{(-\lambda_{\text{B,ref}} T_{\text{ref}} - \lambda_{\text{B,est}} T_{\text{test}})} \\
& \times f_{\text{kibble}}(\lambda_{\text{B,ref}}, \lambda_{\text{B,est}} | \mathbf{a}_B, \mathbf{b}_B, \rho) \\
& \times \left[ \prod_{k=0}^{K-1} (\lambda_{\text{B,ref}} \mathbf{x}_B[k])^{m_{\text{ref}}[k]} (\lambda_{\text{B,est}} \mathbf{x}_B[k])^{m_{\text{test}}[k]} \right] \\
& \times \text{FPT}(\mathbf{x}_B | \mathbf{S}_B)
\end{split}
\tag{21}
\end{equation}
$$

with $M_{\text{ref}} = \sum_{k=0}^{K-1} m_{\text{ref}}[k]$ and $M_{\text{test}} = \sum_{k=0}^{K-1} m_{\text{test}}[k]$.

The priors for vector **x**<sub>B</sub> and the scalar doublet (λ<sub>B,ref</sub>, λ<sub>B,est</sub>) are mutually independent, so the *a priori* probability P(**m**<sub>ref</sub>, **m**<sub>test</sub> | H<sub>0</sub>) may be obtained by first marginalizing over Eq. (21) according to **x**<sub>B</sub>, followed by a second marginalization according to λ<sub>B,ref</sub> and λ<sub>B,est</sub>.

Marginalizing the joint likelihood in Eq. (21) according to **x**<sub>B</sub> requires explicating FPT(**x**<sub>B</sub> | **S**<sub>B</sub>) with the beta distribution B and the dyadic cutting coefficients (b<sub>l,k</sub>)<sub>0≤l≤L,0≤k≤K−1</sub>. Moreover, the derivation requires introducing decomposition coefficients (h<sub>l,k</sub>)<sub>0≤l≤L,0≤k≤K−1</sub> for the couting histogram formed by the sum of **m**<sub>ref</sub> and **m**<sub>test</sub>. These coefficients are calculated as follows:

$\forall k \in [0, K-1], h_{l,k} = m_{ref}[k] + m_{test}[k]$. Then, recursively, from bottom to top of a tree of counts whose structure is given by FPT(**x**<sub>B</sub> | **S**<sub>B</sub>):

$$
\forall l \in [1, L-1], \forall k \in [0, K_{l-1}], h_{l-1,k} = h_{l,2k} + h_{l,2k+1} \quad (22)
$$

After completing all computations, the marginalizations according to **x**<sub>B</sub>, λ<sub>B,ref</sub>, and λ<sub>B,est</sub> yield the following expression for the joint conditional probability of the reference and tested couting spectra under H<sub>0</sub>:
$$ \begin{aligned} \mathbb{P}(\mathbf{m}_{\text{ref}}, \mathbf{m}_{\text{test}} | H_0) = & \frac{T_{\text{ref}}^{M_{\text{ref}}} T_{\text{test}}^{M_{\text{test}}}}{\prod_{k=1}^{K} m_{\text{ref}}[k]! m_{\text{test}}[k]!} \left( \frac{b_B^{a_B}}{\Gamma(a_B)} \right)^2 \frac{\Gamma(a_B + M_{\text{ref}}) \Gamma(a_B + M_{\text{test}})}{(1 - \rho)^{a_B} \tilde{T}_{\text{ref}}^{a_B + M_{\text{ref}}} \tilde{T}_{\text{test}}^{a_B + M_{\text{test}}}} \\ & \times (1 - z)^{-(a_B + M_{\text{ref}} + M_{\text{test}})} {}_2F_1(-M_{\text{ref}}, -M_{\text{test}}, a_B, z) \times \prod_{l=1}^{L} \prod_{k=0}^{K_{l-1}-1} \frac{B(h_{l,2k} + b_{l,2k}, h_{l,2k+1} + b_{l,2k+1})}{B(b_{l,2k}, b_{l,2k+1})} \end{aligned} \quad (23) $$

where parameters $a_B$, $\rho$, $T_{\text{ref}}$, $T_{\text{test}}$, $\tilde{T}_{\text{ref}}$, $\tilde{T}_{\text{test}}$, and $z$ have the same meanings as defined in Section 3.1.

In a similar fashion, we express the joint likelihood of $\mathbf{m}_{\text{ref}}$, $\mathbf{m}_{\text{test}}$, $\mathbf{x}_B$, $\lambda_{B_{\text{ref}}}$, $\lambda_{B_{\text{test}}}$, $\mathbf{x}_S$, and $\lambda_S$ under the alternative hypothesis $H_1$ as:

$$ C_{L-1,k,j} = \sum_{\substack{g+d=j \\ g \le m_{\text{test}}[2k] \\ d \le m_{\text{test}}[2k+1] \\ + b_{L,2k}, h_{L,2k+1} - d + b_{L,2k+1}}} C_{L,2k,g} C_{L,2k+1,d} B(g+s_{L,2k}; d+s_{L,2k+1}) B(h_{2k}-g + b_{L,2k}, h_{L,2k+1}) $$

$$ \begin{aligned} f(\mathbf{m}_{\text{test}}, \mathbf{m}_{\text{ref}}, \lambda_S, \mathbf{x}_S, \lambda_{B_{\text{ref}}}, \lambda_{B_{\text{test}}}, \mathbf{x}_B | H_1) = & \frac{T_{\text{ref}}^{M_{\text{ref}}} T_{\text{test}}^{M_{\text{test}}}}{\prod_{k=1}^{K} m_{\text{ref}}[k]! m_{\text{test}}[k]!} e^{-\lambda_{B_{\text{ref}}} T_{\text{ref}} - (\lambda_S + \lambda_{B_{\text{test}}}) T_{\text{test}}} \\ & \times \frac{b_S^{a_S}}{\Gamma(a_S)} \lambda_S^{a_S-1} e^{-b_S \lambda_S} \times f_{\text{Kibble}}(\lambda_{B_{\text{ref}}}, \lambda_{B_{\text{test}}}|a_B, b_B, \rho) \\ & \times \left[ \prod_{k=0}^{K-1} (\lambda_{B_{\text{ref}}} \mathbf{x}_B[k])^{m_{\text{ref}}[k]} (\lambda_S \mathbf{x}_S[k] + \lambda_{B_{\text{test}}} \mathbf{x}_B[k])^{m_{\text{test}}[k]} \right] \times \text{FPT}(\mathbf{x}_S | \mathbf{S}_S) \times \text{FPT}(\mathbf{x}_B | \mathbf{S}_B) \end{aligned} \quad (24) $$

As before, the marginalizations according to vectors $\mathbf{x}_B$ and $\mathbf{x}_S$ on the one hand, and scalars $\lambda_{B_{\text{ref}}}$, $\lambda_{B_{\text{test}}}$, and $\lambda_S$ on the other hand, are carried out separately.

Marginalizing the joint likelihood in Eq. (24) according to $\mathbf{x}_B$ and $\mathbf{x}_S$ is done by explicitly defining the intermediate simplexes $(\mathbf{x}_{B_{n,k}})_{1 \le n \le L, 0 \le k \le K-1}$ and $(\mathbf{x}_{S_{n,k}})_{1 \le n \le L, 0 \le k \le K-1}$ associated with every level of $\text{FPT}(\mathbf{x}_B | \mathbf{S}_B)$ and $\text{FPT}(\mathbf{x}_S | \mathbf{S}_S)$ respectively:

- Until obtaining the top coefficients $(C_{0,0,j})_{0 \le j \le m_{\text{test}}[k]}$, which we rename as $(\zeta_j)_{0 \le j \le m_{\text{test}}[k]}$ for convenience.

After completing all calculations, the marginalizations with respect to $\mathbf{x}_B$, $\mathbf{x}_S$, $\lambda_{B_{\text{ref}}}$, $\lambda_{B_{\text{test}}}$, and $\lambda_S$ result in the following expression for the joint conditional probability of the reference and tested spectra under $H_1$:

$$ \begin{aligned} \mathbb{P}(\mathbf{m}_{\text{ref}}, \mathbf{m}_{\text{test}} | H_1) = & \frac{T_{\text{ref}}^{M_{\text{ref}}} T_{\text{test}}^{M_{\text{test}}}}{\prod_{k=1}^{K} m_{\text{ref}}[k]! m_{\text{test}}[k]!} \left[ \prod_{l=1}^{L} \prod_{k=0}^{K_{l-1}-1} \frac{1}{B(s_{l,2k}, s_{l,2k+1})} \frac{1}{B(b_{l,2k}, b_{l,2k+1})} \right] \frac{1}{\Gamma(a_S)} \\ & \times \left( \frac{b_S}{T_{\text{test}} + b_S} \right)^{a_S} \times \left( \frac{b_B^{a_B}}{\Gamma(a_B)} \right)^2 \frac{\Gamma(a_B + M_{\text{ref}})}{(1 - \rho)^{a_B} \tilde{T}_{\text{ref}}^{a_B + M_{\text{ref}}} \tilde{T}_{\text{test}}^{a_B + M_{\text{test}}}} \\ & \times (1 - z)^{-(a_B + M_{\text{ref}} + M_{\text{test}})} \times \sum_{j=0}^{M_{\text{test}}} \zeta_j \Gamma(j + a_S) \Gamma(M_{\text{test}} - j + a_B) {}_2F_1(-M_{\text{ref}}, -M_{\text{test}} + j, a_B, z) \psi^j \end{aligned} \quad (26) $$

$$ \begin{aligned} \forall n \in [1, L], \quad \forall k \in [0, K-1], \quad & \mathbf{x}_{B_{n,k}} = \prod_{l=1}^{n} p_{B_{L,\Delta(n-l,k)}} \\ & \mathbf{x}_{S_{n,k}} = \prod_{l=1}^{n} p_{S_{L,\Delta(n-l,k)}} \end{aligned} \quad (25) $$

where $p_{B_{L,\Delta(n-l,k)}}$ and $p_{S_{L,\Delta(n-l,k)}}$ carry the same significance as the introduced in Section 3.2.

After sorting and integrating the expression resulting from the introduction of these simplexes, we observe that marginalization with respect to $\mathbf{x}_B$ and $\mathbf{x}_S$ can be achieved by identifying a series of polynomial coefficients denoted as $(C_{L,k,j})_{1 \le n \le L, 0 \le k \le K_{l-1}, 0 \le j \le m_{\text{test}}[k]}$. These coefficients are associated with the intermediate simplexes to be summed. The procedure for obtaining these coefficients is outlined as follows:

- First, we define the bottom coefficients

$$ \forall k \in [0, K-1], C_{L,k,j} = \begin{pmatrix} m_{\text{test}}[k] \\ j \end{pmatrix}; $$

- Then, we recursively calculate the coefficients for the upper-stage, starting with:

$$ \forall k \in [0, K-1], \forall j \in [0, K_{l-1}], $$

where $\psi$ maintains the same definition as outlined in Section 3.1.

The posterior probabilities of hypotheses $H_0$ and $H_1$ given the reference and tested spectra are obtained determined using Bayes' theorem, following the methodology described in Section 3.1. The Bayes factor for the BKPT test, which plays the same role in the formulation of the statistical test as in aforementioned section, is labeled $\kappa_{\text{BKPT}}$ and is made explicit using Eqs. (23) and (26) as:

$$ \begin{aligned} \kappa_{\text{BKPT}} = & \left[ \prod_{l=1}^{L} \prod_{k=0}^{K_{l-1}-1} \frac{1}{B(s_{l,2k}, s_{l,2k+1})} B(h_{l,2k} + b_{l,2k}, h_{l,2k+1} + b_{l,2k+1}) \right] \\ & \times \left( \frac{b_S}{T_{\text{test}} + b_S} \right)^{a_S} \frac{1}{\Gamma(a_S) \Gamma(a_B + M_{\text{test}}) {}_2F_1(-M_{\text{ref}}, -M_{\text{test}}, a_B, z)} \\ & \times \sum_{j=0}^{M_{\text{test}}} \zeta_j \Gamma(j + a_S) \Gamma(M_{\text{test}} - j + a_B) {}_2F_1(-M_{\text{ref}}, -M_{\text{test}} + j, a_B, z) \psi^j \end{aligned} \quad (27) $$

Eventually, the BKPT hypothesis test against a given $\alpha$-risk (or FAR) is formulated as follows:

- If $\frac{1}{1+\kappa_{\text{BKPT}}} < \eta_\alpha$, then $H_0$ is rejected and a uranium contamination signal is detected;

- Otherwise, $H_0$ is accepted.
Again, the relation between the decision threshold $\eta_{\alpha}$ and the $\alpha$ - risk must be determined numerically using a Monte Carlo method.

## Parametrization of the FPT model

The gamma distribution of Kibble is parametrized as described in Section 3.1 for the BKD test. Following the procedure outlined in Section 3.2.1, we ensure that the means of $FPT(\mathbf{x}_B|\mathbf{S}_B)$ and $FPT(\mathbf{x}_S|\mathbf{S}_S)$ match the *a priori* signatures of the radiological background and a uranium contamination signal with a known U-235 enrichment level, respectively. Consequently, the scale factors $c_0$ and $\sigma$ are the only parameters remaining to be selected. It is advisable to choose a value close to 1 for $c_0$ to limit the influence of the prior on the data at the first level of the Pólya tree.

Recalling that for every channel $k$ in the spectral ROI, $b_{L,k} = s_B[k] \bullet (2\sigma)^L$ and $s_{L,k} = s_S[k] \bullet (2\sigma)^L$, we observe that increasing $\sigma$ results in a stricter adherence of the FPT to the *a priori* signatures $s_B$ and $s_S$, thereby making the models more rigid. A satisfactory tradeoff between the information provided by these signatures and the robustness of the test against variability in the background and actual U-235 enrichment level was numerically determined to be achieved at $\sigma = 2$.

To conclude this section, it is noteworthy that the FPT model for $\mathbf{x}_B$ and $\mathbf{x}_S$ converges toward the Dirichlet model used in the BKD test when $\sigma = \frac{1}{2}$ and $c_0 = K \cdot \nu_S$, as per the notations in Section 3.1. Generally, for a Dirichlet prior, it holds that $K \cdot \nu_S \gg 1$.

## Results and discussion

The numerical results presented in this section use ROC curves to compare the performance of the BKPT against various hypothesis tests. These tests include frequentist tests such as the KS, the cumulative function of a negative binomial law (CNB), and the cumulative function of a beta binomial law (CBB). Additionally, Bayesian tests like the Bayesian absolute test by mixture of multivariate Poisson laws (BAM), and Bayesian relative test by mixture of multinomial laws (BRM), and BKD are also considered. Random tested spectra were generated and ROC curves constructed using a Python-based platform. As outlined in Section 3.2.2, implementing the BKPT test requires setting the parameters $c_0$ and $\sigma$, which influence the FPT model towards either a smooth or Dirichlet-like prior. The results for these two parameters are detailed in the following sections, under the specified configurations:

- Bayesian hypothesis tests: BAM, BRM, BKD and BKPT;

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** ROC curves for an enrichment level of Nat wt%, $\sigma = 1/2$, $c_0 = 10$ in BKPT, an integration time $t = 1000$ s, and (a) stationary, or (b) unstationary radiological backgrounds.</figcaption></figure>

- Frequentist hypothesis tests: KS, CNB and CBB;

- Uranium enrichment levels: Nat wt%, 8 wt%;

- Examined subunits: SU1 and SU3;

- Source: Uranium gamma-ray peaks contamination;

- Detector: HPGe detector;

- Source-detector distance: $d = 22.36$ cm for a uranium surface activity of 500 Bq.0.1 m⁻² over $2\pi$;

- BKPT parameters: $\sigma = \{2, 1/2\}$ and $c_0 = \{1, 10\}$;

- BKD parameter: $\xi = 1$;

- Weighting factors: $\nu_S = \nu_B = 10$;

- Integration time: $t = \{1000, 3000, 5000\}$ s;

- Number of iterations $N = 1000$;

- Radiological background: stationary and unstationary with respect to shape and intensity.

It is important to highlight that the implementation of the BKD test with respect to the Pearson correlation coefficient $\rho$ follows the same approach as outlined in [9]. The influence of $\rho$ on the BKD test was calibrated using an intermediate parameter $\xi$, a parametrization justified in article [9]. In this context, we chose $\xi = 1$, as it provides an intermediate solution between the absolute (BAM) and relative (BRM) tests.

Previously, we mentioned that the normalized prior vectors $s_B$ and $s_S$, fundamental to all three Bayesian tests discussed, were constructed using in-situ data acquisitions of background spectra and experimentally-calibrated MCNP6.1 simulations with Eu-152 and variously-enriched U samples. This approach ensures their empirical relevance. Considering the context of a surface with relatively low U-235 and U-238 activity, our proof-of-concept study requires rigorous control over the isotopic signal during detection trials. This signal must be compared with MSAD, and numerous trials must be conducted with repeatability assurance at this stage. Therefore, we opted for a detailed numerical study rather than a large-scale field campaign. However, a measurement campaign is scheduled to follow, involving the acquisition of tested spectra from concrete surfaces suspected of containing uranium traces, subsequent to the presentation of the results in this section.

## Performance study of the BKPT test for a dirichlet-like FPT model

Throughout Section 3.2, we introduced the BKPT detection test as a generalization of the BKD test, which should be able to replicate the features of the latter with appropriate parameters, namely $\sigma = 1/2$ and
<figure><img src="image_7.png" /><figcaption>Fig. 5. ROC curves for an enrichment level of 8 wt%, with $\sigma = 1/2$ and $c_0 = 10$ in BKPT, and an integration time of $t = 1000$ s, for (a) stationary and (b) unstationary radiological backgrounds.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. ROC curves for a natural (Nat wt%) enrichment level, with $\sigma = 2$ and $c_0 = 1$ in BKPT, and an integration time of $t = 3000$ s, for (a) stationary and (b) unstationary radiological backgrounds.</figcaption></figure>

$c_0 = K_{U,S,B}$ 10. To verify this claim, we present ROC curves for both the BKD and the BKPT test, parametrized with $\sigma = 1/2$, $c_0 = 10$, for natural (Fig. 4) and enriched U (Fig. 5). These curves are shown against both stationary (Figs. 4a and 5a) and unstationary (Figs. 4b and 5b) radiological backgrounds.

From Figs. 4-5, we observe that, regardless of the radiological background's stationarity status or the U-235 enrichment level, both BKPT and BKD tests exhibit nearly identical detection behavior. This observation explicitly confirms the convergence of the FPT model towards the Dirichlet model, as indicated in Section 3.2.2. This result corroborates that the BKPT test is a generalized version of the BKD test, demonstrating sensitivity to both the shape and density of the spectrum. Consequently, it can be concluded that the BKPT test provides a comprehensive detection process that accommodates the diverse requirements and constraints that may arise in the context of a large decommissioning project.

Performance estimation of the BKPT test for a smooth FPT model

In this section, ROC curves for the BKPT and BKD tests are presented, retaining $\sigma = 2$ and $c_0 = 1$ as parameters for the FPT model implemented in the latter. The curves are plotted for natural (Fig. 6) and enriched U (Figs. 7 and 8), and against both stationary (Figs. 6a, 7a, 8) and unstationary (Figs. 6b and 7b) radiological backgrounds.

Figs. 6-8 illustrate a statistically significant improvement in detection performance when moving from an optimally calibrated ($\xi = 1$) BKD test to a smooth ($\sigma = 2$) FPT model in the BKPT test. As an illustration, under a stationary radiological background and enrichment levels of Nat wt% and 8 wt%, BKPT yielded respective FOMs $\max_{ROC}(\text{TDR} - \text{FAR}) = (0.62 \pm 0.01)$ (extracted from Fig. 6(a)) and $(0.90 \pm 0.01)$ (from Fig. 7(a)), which compare favorably to the respective figures $(0.59 \pm 0.01)$ and $(0.87 \pm 0.01)$ obtained with BKD.

Additionally, as shown in Fig. 8, it is evident that the performance of BKPT closely aligns with that of an ideal detection test, specifically one with a 1% significance level, under conditions characterized by
<figure><img src="image_9.png" /><figcaption>Fig. 7. ROC curves for an enrichment level of 8 wt%, with $\sigma = 2$ and $c_0 = 1$ in BKPT, and an integration time of $t = 3000$ s, for (a) stationary and (b) unstationary radiological backgrounds.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. ROC curves for an enrichment level of 8 wt%, with $\sigma = 2$ and $c_0 = 1$ in BKPT, and an integration time of $t = 5000$ s, for stationary radiological backgrounds.</figcaption></figure>

increased enrichment and prolonged integration time.

Another noteworthy observation apparent in this set of curves is the adverse impact of a slightly unstationary background on the achievable tradeoff between TDR and FAR for BKPT. For instance, there is a noticeable degradation in sensitivity for an enrichment level of 8 wt% at $t = 3000$ s, from $\max_{ROC}(\text{TDR} - \text{FAR}) = (0.90 \pm 0.01)$ in Fig. 7(a) down to TDR = (0.70 ± 0.01)% in Fig. 7(b).

Performance study of BKPT compared to frequentist and Bayesian
hypothesis tests

In this section, ROC curves for the BKPT and previous related work, including frequentist tests (KS, CNB, CBB) and Bayesian tests (BAM, BRM, BKD) are presented. In this scenario, Dirichlet-like a priori ($\sigma = 1/2$, $c_0 = 10$) and smooth ($\sigma = 2$, $c_0 = 1$) parameters are retained for the FPT model.

The reproduced ROC curves provide evidence of the superior per-
formance of the BKPT over the frequentist tests, regardless of the
parameterization of the FPT model. As depicted in Fig. 9(a), for an

enrichment level of 8 wt% and $\sigma = 1/2$, $c_0 = 10$, the BKPT and CNB tests provide FOM $\max_{ROC}(\text{TDR} - \text{FAR})$ values of $(0.69 \pm 0.01)$ and $(0.21 \pm 0.03)$, respectively. As expected, BKPT with $\sigma = 1/2$, $c_0 = 10$ exhibits a similar TDR/FAR behavior to its Bayesian counterpart in Fig. 9(a), while BKPT with $\sigma = 2$, $c_0 = 1$ demonstrates a notable increase in the FOM in Fig. 9(b). For example, BKPT achieves a FOM $\max_{ROC}(\text{TDR} - \text{FAR}) = (0.72 \pm 0.01)$, compared to FOMs of $(0.68 \pm 0.01)$ and $(0.20 \pm 0.03)$ with BAM and CNB, respectively.

These results highlight the flexibility and comprehensiveness of the newly introduced BKPT statistical test for a large-scale surface characterizations within walls. Even with a detector with as low an energy resolution as an HPGe detector, there is notable gain improvement in detection performance by considering interchannel counting correlations. Additionally, there is no clear preference for intensity background variability over shape variability, suggesting potentially significant performance differences between BKD and smooth BKPT tests when applied to spectra from lower-energy resolution detectors like NaI(Tl) scintillators.

However, there are inherent limitations to overcome. Bayesian pro-
cedures do not support integration times below the order of 10³ seconds.
Considering the substantial levels of uranium FEP and background
within the associated ROI, achieving a SNR of about 3.5 at this order
remains highly challenging in nuclear counting applications [21,22].
Further, irregularities on the concrete surface, leading to variable dis-
tances between the source and detector, poses significant challenges in
accurately quantifying uranium activity. To enhance the fieldability of
our approach, integrating Bayesian inference with real-time Convolu-
tional Neural Networks (CNNs) could be explored in the future for real-
world applications.

Furthermore, the parametrization of the FPT model is crucial for
adapting the BKPT approach to field applications. In practical terms, for
a spectrum composed of a Compton baseline and peaks, we recommend
setting $\sigma = 2$. Increasing this value enables us to closely align with $s_S$ (
the normalized a priori spectral signature of the uranium signal), while
decreasing it allows us to deviate further from $s_S$.

Conclusion

The objective of this paper was to develop an effective method based
on Bayesian decision-making for detecting low-level radioactivity
against challenging external measurement conditions, particularly
against a standard radiological background, within the field of radiation
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** ROC curves for an enrichment level of 8 wt%, with stationary radiological backgrounds, and an integration time of $t = 1000$ s, comparing scenarios (a) with $\sigma = 1/2$ and $c_0 = 10$, and (b) $\sigma = 2$ and $c_0 = 1$ in BKPT test.</figcaption></figure>

detection. Detecting a weak uranium signal, characterized by a very low SNR, amid a highly variable radiological background presents significant operational and technical constraints in selecting suitable detection methods. Considering these challenges, an innovative Bayesian detection approach based on a Pólya tree-type *a priori*, termed the BKPT test, designed to address this complex detection problem. The novelty of the BKPT test primarily relies on the FPT, enabling for the modeling of a higher degree of correlations between counts recorded in neighboring spectral channels. The implementation of the BKPT test involves setting parameters $c_0$ and $\sigma$, which guide the FPT model towards either a smooth or Dirichlet-like priori. This adjustment enhances the test, making it more comprehensive within a Bayesian framework.

Based on ROC curves, we have confirmed that the BKPT test achieves comparable FOM to the BKD test when its underlying FPT model assumes a Dirichlet-like parametrization. Conversely, with a smooth FPT parametrization in BKPT, we observed a superior tradeoff between specificity (or FAR) and sensitivity (or TDR) compared to its discrete counterpart BKD. For instance, the optima in associated ROC curves increased by 3% to 28% against unstationary and stationary radiological backgrounds, respectively. Moreover, this improved tradeoff closely approximates, within a 1% margin, the performance of an ideal detection test under conditions of high enrichment level and extended integration time. The superiority of BKPT is further validated through comparisons with previous related works, including absolute (BAM) and relative (BRM) Bayesian hypothesis tests, as well as more significant gains observed in absolute (CNB) and relative (CBB) detection tests. For instance, when considering a stationary radiological background, the estimated gains are 52% and 4% with BAM and CNB, respectively.

The ability of the Bayesian statistical tests to discriminate the spectral signature of a low uranium signal against the variable background of a closed, concrete-based environment compares advantageously to analogous frequentist tests such as KS, CNB and CBB. This emphasizes the merits of Bayesian inference in decision-making process, particularly in scenarios involving radiological alarms. Bayesian methods not only help in raising alarms but also provide crucial supplementary information for decision-makers, enhancing risk management strategies and mitigating environmental impacts. This demonstrates the robustness and applicability of Bayesian approaches in complex detection and monitoring tasks.

The remarkable performance observed will undoubtedly inspire researchers to innovate and develop new statistical methods, leveraging both parametric and non-parametric Bayesian modeling techniques. A particular focus will be on applying the BKPT statistical test to

experimental uranium spectra to enhance the interpretation of simulation-based ROC curves and time-dependent detection rates. This is especially crucial in the demanding context of large-surface decommissioning projects. Collaborative initiatives with the Laboratoire National Henri Becquerel (LNHB, CEA List) are aimed at advancing metrological manufacturing and characterization of uniformly dispersed uranium sources. These efforts are aligned with the development of regulatory MSAD, emphasizing the importance of accurate and reliable detection methods in radiological monitoring and control. Another promising research avenue to explore involves applying the BKPT test to conduct line-by-line analysis of U-235.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## References

[1] Petraglia A, Sirignano C, Buompane R, D'Onofrio A, Esposito AM, Terrasi F, et al. Space-time Bayesian analysis of the environmental impact of a dismissing nuclear power plant. J Environ Radioact 2020;218:106241. https://doi.org/10.1016/j.jenvrad.2020.106241.

[2] Kim J, Lim KT, Ko K, Ko E, Cho G. Radioisotope identification and nonintrusive depth estimation of localized low-level radioactive contaminants using Bayesian inference. Sensors 2020;20:95-108. https://doi.org/10.3390/s20010095.

[3] Guo C, Khan F, Imtiaz S. Copula-based Bayesian network model for process system risk assessment. Process Saf Environ Prot 2019;123:317-26. https://doi.org/10.1016/j.psep.2019.01.022.

[4] Jetti HV, Ferrero A, Salicone S. A modified Bayes' theorem for reliable conformity assessment in industrial metrology. Measurement 2021;184:109967. https://doi.org/10.1016/j.measurement.2021.109967.

[5] Kim J, Lim KT, Kim J, Kim Y, Kim H, Cho G. Quantification and uncertainty analysis of low-resolution gamma-ray spectrometry using Bayesian inference. Nucl Instrum Meth A 2020;953:163144. https://doi.org/10.1016/j.nima.2019.163144.

[6] Michaud IJ, Schmidt K, Smith RC, Mattingly J. A hierarchical Bayesian model for background variation in radiation source localization. Nucl Instrum Meth A 2021; 1002:165288. https://doi.org/10.1016/j.nima.2021.165288.

[7] Hasan MM, Rogiers B, Laloy E, Camps J, Rutten J, M., Huysmans, Estimation of soil radioactivity-depth profiles using Bayesian inversion of borehole gamma spectrometry data. J Environ Radioact 2023;257:107077. https://doi.org/10.1016/j.jenvrad.2022.107077.

[8] Knoll GF. Radiation Detection and Measurement. 4th ed. Hoboken, NJ, USA: Wiley; 2010.

[9] Arahmane H, Dumazert J, Barat E, Dautremer T, Carrel F, Dufour N, et al. Bayesian Inference based on a Bivariate Gamma Distribution of Kibble for Low-Level Radioactivity Detection in nuclear decommissioning operations. Process Saf Environ Prot 2022;163:727-42. https://doi.org/10.1016/j.psep.2022.05.034.

[10] Arahmane H, Dumazert J, Barat E, Dautremer T, Dufour N, Carrel F, et al. An absolute and relative Bayesian method for nuclear decommissioning: low-level
radioactivity measurement with gamma-ray spectrometry. IEEE Trans Instrum
Meas 2021;70. https://doi.org/10.1109/TIM.2021.3065419.

[11] Dufour N, Dumazert J, Barat E, Bertrand GHV, Carrel F, Dautremer T, et al. Measurement of Low-Activity Uranium Contamination by Gamma-Ray Spectrometry for Nuclear Decommissioning. Nucl Instrum Meth A 2020;951. https://doi.org/10.1016/j.nima.2019.162976.

[12] D.B. Pelowitz MCNP6TM User's Manual Version 1.0 2013 Los Alamos National Laboratory, USA CP-13-00364, Tech. Rep. LA-807.

[13] Metz CE. Basic principles of ROC analysis. Seminars Nucl Med 1978;8:283-98.
https://doi.org/10.1016/s0001-2998 (78)80014-2.

[14] Bé MM, Chisté V. Table de radionucléides, LNE-LNHB/CEA, Paris-Saclay. France: Tech. Rep; 2007.

[15] GMX Series Coaxial HPGe Detector Product Configuration Guide.[Online].
Available: https://www.ortec-online.com/-/media/ametekortec/brochures/
gamma-x.pdf.

[16] Tsutsumi M, Oishi T, Kinouchi N, Sakamoto R, Yoshida M. Simulation of the background for gamma detection system in the indoor environment of concrete

buidings. J Nucl Sci Technol 2001;38:1109-14. https://doi.org/10.1080/
18811248.2001.9715143.

[17] Balakrishna N, Lai C-D. Continuous Bivariate Distributions. second ed. New York: Springer; 2009.

[18] Chen L-S, Tzeng I-S, Lin C-T. Bivariate generalized gamma distributions of kibble's type. Statistics 2014;48:933-49. https://doi.org/10.1080/02331888.2012.760092.

[19] Lavine M. Some aspects of Pólya tree distributions for statistical modelling. The Annals of Statistics 1992;20:1222-35. https://doi.org/10.1214/ aos/1176348767.

[20] Lavine M. More aspects of polya tree distributions for statistical modelling. The Annals of Statistics 1994;22:1161-76. https://doi.org/10.1214/ aos/1176325623.

[21] Dumazert J, Coulon R, Kondrasovs V, Boudergui K, Moline Y, Sannié G, et al. A robust hypothesis test for the sensitive detection of constant speed radiation moving sources. Nucl Instrum Meth A 2015;795:335-42. https://doi.org/10.1016/j.nima.2015.06.016.

[22] Arahmane H, Dumazert J, Barat E, Dautremer T, Carrel F, Dufour N, et al. Statistical approach for radioactivity detection: A brief review. J Environ Radioact 2024;272. https://doi.org/10.1016/j.jenvrad.2023.107358.