## 3 Numerical Results  
### Dataset: Measured Bottarga Samples  
It has been already mentioned that the bottarga samples come from four different areas: two of them located in Greece namely Messolongi and Preveza. Both areas are producing Bottarga products, however only the Bottarga originating in Messolongi has been registered since 1996 as a Protected Designations of Origin product (PDO-GR-0446) with the name \"Avgotaracho Messolongiou\". The rest of samples come from two other regions of the world (Australia and Mauritania). The samples were labeled by codes \"A1\" to \"A4\", see Table 1. Bottarga production in the geographical region of Messolongi is very limited and consequently expensive; thus, a random sampling scheme of five samples (derived from batches with different lot numbers) has been used. For regions outside Messolongi, a random selection performed on a retailer level trying to pick up different samples coming from different lot number batches. The measured mean values of each area are given in the following Table 2. The scatter plots of the mean value of each sample are shown in Figure 1 (each dot corresponds to the mean values of a sample), while the main diagonal is the kernel density estimation of each pair of the primary features. A clear separation of bottarga samples (orange dots) from Messolongi is observed. Table 3 shows the pairwise Pearson correlation coefficients of the mean values of all samples [41]. It is obvious that there is a very strong linear association between all combinations of primary features; this is a strong indication of multicollinearity [34]. A practical interpretation would be that additional ratios of stable isotopes might not add predictive power to the model. By using the family of models described in Equations (1)-(9), the data (i.e., mean value of subsamples per sample), and MCMC/HMC (software framework PYMC3 [36]), we estimate the posterior probability distributions of all candidate models in Appendix A. The sampler of choice is NUTS [38], we use two chains for each model, and we draw 15,000 parameter samples (after 3000 warming up samples). Based on the model comparison methodology of Section 2.4, the proposed probabilistic machine learning model for the characterization of protected designation of origin of bottarga from Messolongi is the following one:  
\[m_{\rm 534 SV-CDT\ (\%_{\omega})^{2}}\sim\rm N(0,\infty) \tag{13}\]  
\[\rm c\ \sim\rm N(0,\infty) \tag{14}\]  
\[y_{i}\sim B(n,\frac{1}{1+e^{-m_{\rm 534 SV-CDT\ (\%_{\omega})^{2}}\approx 34 SV-CDT\ (\%_{ \omega})^{2}+c}}) \tag{15}\] \begin{table}
\begin{tabular}{l c c c c} \hline \hline
**Code** & **Area** & **Mean \(\rm\xi^{15}N_{AIR}\) (\%\(\%\))** & **Mean \(\rm\xi^{13}C_{V\_TDB}\)** & **Mean \(\rm\xi^{34}S_{V\_CDT}\) (\%\(\%\))** \\ \hline A1 & Mauritania & 9.6442 (0.0073) & \(-19.0075\) (0.0084) & 16.0344 (0.0109) \\ A2 & Preveza Greece & 10.8551 (0.0146) & \(-18.0034\) (0.0187) & 15.1658 (0.0073) \\ A3 & Messolongi Greece & 7.5246 (0.0056) & \(-18.4259\) (0.0068) & 11.9966 (0.0135) \\ A4 & Australia & 5.6429 (0.0050) & \(-14.5148\) (0.0088) & 4.5341 (0.0094) \\ \hline \hline \end{tabular}
\end{table} 
Table 2: Bottarga Analysis Results (mean \(\delta\) values and standard error). \begin{table}
\begin{tabular}{c c c c} \hline \hline  & \(\rm\delta^{15}N_{AIR}\) (\%\(\%\)) & \(\rm\delta^{13}C_{V\_TDB}\) (\%\(\%\)) & \(\rm\delta^{34}S_{V\_CDT}\) (\%\(\%\)) \\ \hline \(\rm\delta^{15}N_{AIR}\) (\%) & 1.00 & \(-0.85\) & 0.95 \\ \(\rm\delta^{13}C_{V\_TDB}\) (\%) & \(-0.85\) & 1.00 & \(-0.97\) \\ \(\rm\delta^{34}S_{V\_CDT}\) (\%) & 0.95 & \(-0.97\) & 1.00 \\ \hline \hline \end{tabular}
\end{table} 
Table 3: Pearson coefficient correlation matrix of all measured bottarga samples. This is a model that uses the squared value of \(\delta^{34}\)S\({}_{\text{V-CDT}}\) (%) and an intercept in the linear part of the logistic function (15) to define the probability of success of the Bernoulli trials.