Inferring the shape of data: A probabilistic framework for analyzing experiments in the natural sciences

To begin, we consider the problem of specifying a model of the latent structure of a dataset for the purpose of mimicking visual recognition. In our framework, a dataset, $\bm{y}$, is a tensor whose components are the individual, scalar datapoints. Regardless of the experimental relationships between those datapoints (i.e., the organisational structure of the tensor), for simplicity we can reshape $\bm{y}$ into an N-dimensional vector $\bm{y}=[y_{1},...,y_{N}]$, where $y_{i}$ are the scalar datapoints. One can imagine a dataset $\bm{y}$ collected using a particular instrument, in a particular location, on a particular day, and by a particular researcher. Altogether, these specific factors might induce systematic differences in $\bm{y}$ relative to an otherwise equivalent experiment. For example, an optical filter in an instrument might slowly oxidise, which could reduce the intensity of light incident upon the detector and, over a period of months, yield $\bm{y}$ with different scales (i.e., units). Similarly, overhead lights might be left on accidentally when making an optical measurement, which could increase the background photons incident on the detector and yield $\bm{y}$ with different relative offsets. Likewise, local vibrations might vary from day to day, which could affect the stability of an instrument and yield $\bm{y}$ with different amounts of noise. Yet factors such as the scaling, offset, or amount of noise in a measurement generally do not alter the underlying physical processes that give rise to the measured data, and thus should not affect the latent structure (i.e., the shape) of $\bm{y}$. Instead, these factors often act as irregularities, or nuisances, that can limit our ability to model the shape of $\bm{y}$, hence our use of the term ‘nuisance parameters’ to describe them. With this in mind, we define a template, $\bm{x}=[x_{1},...x_{N}]$, as a particular N-dimensional vector of scalar quantities which is related to $\bm{y}$ through the following transformation

In this equation, $m$ and $b$ are nuisance parameters representing changes to scale and offset, respectively, $\bm{\xi}=[\xi_{1},...,\xi_{N}]$ is a nuisance parameter composed of stochastic terms representing the experimental noise, and $N$ is the number of components in $\bm{x}$ or $\bm{y}$. If we recall our definition of shape as correlations within data which derive from fundamental physical processes, we can conceptually understand a template, $\bm{x}$, as an ideal representation of these correlations, without noise or background. To avoid confusion, we note that our definition of shape is distinct from those that take shape to mean a boundary or segmentation in data Dryden and Mardia (2016), and that it is this choice of definition which enables our framework and aligns it with the intuitive visual analyses performed manually by researchers.
It is important to note that the shape of $\bm{y}$, regardless of any distortions caused by the experimental nuisance parameters described above, may often be reasonably described by many different $\bm{x}$s. Indeed, there are no restrictions on what specific $\bm{x}$s one may choose as templates. Different $\bm{x}$s might depend upon different levels of theory, the particular details of the experimental setup, and even sample-to-sample variability. For example, the laws of diffraction dictate the point spread function (PSF) that describes the shape of point emissions in a fluorescence microscopy image Pawley (2006). However, for a standard microscope, the PSF can be modeled by an Airy disk; a two-dimensional circular Gaussian function; or, in order to incorporate an astigmatism correction, even a two-dimensional elliptical Gaussian function Shen et al. (2017). Each of these models of the PSF provides a distinct, theoretically valid $\bm{x}$ capable of modeling the shape of $\bm{y}$ with varying degrees of complexity. Alternatively, one’s $\bm{x}$s could be empirically derived from data previously recorded in other experiments. In the context of the above example, an $\bm{x}$ can be created from fluorescence microscopy images of point emitter-like samples without needing to explicitly invoke a theory of diffraction, and such empirically derived $\bm{x}$s might even model the latent structure of $\bm{y}$ more effectively than analytical, theory-derived $\bm{x}$s. Regardless of the complexity of $\bm{x}$ or its origin, once formulated, it is directly related to a $\bm{y}$ by the simple affine transformation given in Eqn. (1).
The choice of which $\bm{x}$ (or set of $\bm{x}$s) one uses to model the shape of $\bm{y}$ depends not only on the experimental technique but also on the level of precision required for that particular analysis. When using this mathematical framework to analyse experimental data from the natural sciences, one can invoke prior knowledge of the physics governing the experiment which gave rise to the data to constrain the choice of templates used in the analysis. Thus, while templates with higher complexities (e.g., an Airy disk as a model for a PSF) may be required for certain applications, in other cases, a less complex template (e.g., a 2D Gaussian as model for a PSF) can perform just as effectively while greatly reducing the computational cost of the analysis. The flexibility in choice of templates enabled by our framework can greatly increase the efficiency and effectiveness of an analysis method (see Sec. III), however, determining which of the chosen $\bm{x}$s, if any, is the optimal template requires that we first compute how well the shape of $\bm{y}$ is explained by a given $\bm{x}$.

After defining an $\bm{x}$, we quantify the degree to which it describes the latent structure of $\bm{y}$, regardless of the nuisance parameters described above in Eqn. (1). For the $k$^th template, $\bm{x}_{k}$, in a set of templates, $\{\bm{x}\}=\{\bm{x}_{1},...,\bm{x}_{K}\}$, this means calculating a marginalised likelihood probability called the evidence, $P(\bm{y}|\bm{x}_{k},M_{0})$. Here, the conditional $M_{0}$ represents all of the details of the experiment, previous knowledge about the system, and particulars of the analysis method(s)—including which templates have been incorporated into the chosen $\{\bm{x}\}$. The expression for the evidence of $\bm{x}_{k}$ is the marginalisation of the joint probability, which is given by

In Eqn. (2), $p(\bm{y}|\bm{x}_{k},\bm{\xi},m,b,M_{0})$ is called the likelihood, and it represents the probability density of observing $\bm{y}$ for a given $\bm{x}_{k}$ and given values of the nuisance parameters; $p(\bm{\xi},m,b\,|\,M_{0})$ is called the prior, and it represents the joint probability density of those particular nuisance parameter values based on the prior knowledge specified by $M_{0}$.
In this work, we have used combinations of different likelihoods and priors to derive a set of evidences, expressed in closed-form, that are particularly useful for calculating the shape of data in a variety of experimental situations. For all of the cases presented here, we have assumed in our $M_{0}$ that the $\xi_{i}$ are uncorrelated, such that $\langle\xi_{i}\rangle$ = 0, and $\langle\xi_{i},\xi_{j}\rangle$ = $\tau^{-1}\delta_{ij}$, where $\tau$ is a constant called the precision and $\delta_{ij}$ is the Kronecker delta. While this assumption is not a requirement of our approach, this noise model is often experimentally reasonable, and it has allowed us to present analytical solutions to evidence integrals in many general situations (see Supplemental Materials, Section 2 for other noise models). Together with Eqn. (1), this assumption yields the following likelihood function:

Very similar likelihood functions arise with this noise model when $m$ is known to be 0 or 1, and/or $b$ is 0.
Specifying the probability expression for the prior–the second term in the integrand in Eqn. (2)–requires that we mathematically represent our previous knowledge of how $m$, $b$, and $\tau$ are distributed in the experiments of interest Jaynes (2003). In particular, the prior dictates the integration bounds of Eqn. (2) by determining the values that are possible for these parameters to assume (i.e., regions where the prior probability is non-zero). For the results derived here, we have used so-called ‘maximum entropy’ priors, which allow us to encode information and constraints into our prior probability expressions, without dictating their functional form in an ad hoc manner Jaynes (2003). If we assume in $M_{0}$ that we only know that $m$, $b$, and $\tau$ are within some range and that we do not know the magnitude of $\tau$ (i.e., the amount of noise we expect), then the corresponding maximum entropy priors are a uniform distribution for $m$ and $b$, and a uniform distribution over the logarithm of $\tau$ (see Supplemental Materials, Section 1). If we further assume that $m$, $b$, and $\tau$ are independent, then the corresponding joint prior of these parameters within the given ranges is

where the shorthand $\Delta f\equiv f_{max}-f_{min}$ defines the range of a parameter.
In order to analytically integrate Eqn. (2), the integration bounds in the prior must be explicitly defined. We note that a positive transformation of an $\bm{x}_{k}$ ($m>0$) can have a distinct physical interpretation from a negative transformation ($m<0$). Thus, in order to differentiate between these two cases and properly model the underlying shape of $\bm{y}$, we impose that $\bm{x}_{k}$ and $\bm{y}$ be oriented in the same direction. This constraint can be encoded into the calculation by only considering orientation-preserving (i.e., positive scaling) affine transforms of the $\bm{x}_{k}$ in Eqn. (1). To explicitly include this information in the prior, and thus in our $M_{0}$, we therefore use $m_{min}=0$ rather than some $m_{min}<0$. In the case that the negative transformation ($m<0$) is of interest, we note that $-\bm{x}_{k}$ with $m>0$ is equivalent to $\bm{x}_{k}$ with $m<0$. Closed-form expressions for the evidence derived using other integration bounds are also provided in the Supplemental Materials. Additionally, to keep the prior normalised and avoid using a so-called ‘improper’ prior, the minimum and maximum values must be chosen such that $\Delta m$, $\Delta b$, and $\Delta\ln\tau$ are not infinite. For the purposes of a tractable integration Gradshteyn et al. (2014); Ng and Geller (1969), we have used such large negative and positive values that the integration bounds in Eqn. (2) can be approximated as $m\in[0,\infty)$, $b\in(-\infty,\infty)$, and $\tau\in[0,\infty)$. While the exact values of the bounds are important and should be chosen judiciously, we note that the resulting prior normalisation terms end up canceling in subsequent steps during BMS (see below). Using the integration bounds discussed above, the closed-form probability expression for the evidence calculated using Eqn. (2) is

We compare the performance of different templates by using BMS Jaynes (2003); Bishop (2006); Kinz-Thompson et al. (2021) in order to calculate the probability that each $\bm{x}_{k}$ is the best description of the shape of the data, $\bm{y}$. This calculation is conditionally dependent on the assumptions in $M_{0}$, which define the specifics of the analysis method, including the composition of $\{\bm{x}\}$. Multiple distinct analysis methods can consequently be developed by using different $M_{0}$s to tailor their effectiveness to individual experimental situations and systems. For any chosen $M_{0}$, an appropriate template prior probability for $\bm{x}_{k}$, $P(\bm{x}_{k}|M_{0})$, must then be assigned, for example, by: (i) using an equal a priori assignment of $K^{-1}$, where $K$ represents the number of templates in $\{\bm{x}\}$; (ii) learning prior values from separate experiments; or even (iii) using a Dirichlet process or hierarchical Dirichlet process Teh et al. (2006) for a non-parametric ‘infinite’ set of templates. Once all of the $P(\bm{x}_{k}|M_{0})$ have been assigned, Bayes’ rule can be used to perform BMS and compute the template posterior probability as

This expression represents the probability of an $\bm{x}_{k}$ given the observed data $\bm{y}$ and, thus, may be used to identify the $\bm{x}_{k}$ in $\{\bm{x}\}$ that most optimally describes the latent structure of $\bm{y}$ (for a specific choice of $M_{0}$). Using Eqn. (7) is therefore a quantitative means by which the underlying shape of experimental data may be determined. Furthermore, by considering a ‘background’-shaped $\bm{x}_{k}$ and/or just the presence of noise (i.e., Eqn. (6)) in the BMS process, this approach can also validate whether using the most probable $\bm{x}_{k}$ to describe the shape of $\bm{y}$ is justified, or whether the shape of $\bm{y}$ can be better explained as just noise in the data. Altogether, this BMS process sets up an objective, quantitative, researcher-independent metric for not only determining the shape of experimental data, but also validating such shape assignments.
The shape-calculation equations we report above describe a relationship between ideal distributions (i.e., $\bm{x}_{k}$) and noisy signals (i.e., $\bm{y}$) that is independent of many experimental details which would otherwise complicate the analysis being performed. The only requirements are that both $\bm{x}_{k}$ and $\bm{y}$ exist in the same data-space and are vectors of the same size. Practically, however, most templates are generated from some underlying model that exists in a separate ‘model-space’ distinct from the data-space of $\bm{x}_{k}$ and $\bm{y}$. Relating such a model-space to data-space requires that a set of parameters, $\{\theta\}$, be used to map the model to an $\bm{x}_{k}$. For example, a model of a three-dimensional object being projected onto a two-dimensional image may use the Euler angles of the object to generate $\bm{x}_{k}$s with different orientations in the two-dimensional image data-space.
Generally, when using such a model to generate $\bm{x}_{k}$s for identifying the shape of $\bm{y}$, the template posterior probabilities for an entire group of model-associated $\bm{x}_{k}$s must be calculated to account for the many possible ways that the single model could have been mapped into data-space. Having performed all of those calculations, it is then possible to marginalise out the dependence upon some of the $\{\theta\}$ from the model. In the example above, marginalising out the Euler angles would yield the posterior probability that the shape of $\bm{y}$ corresponds to a two-dimensional image of the model, regardless of not knowing the true orientation of the three-dimensional object being projected into the image. Thus, this type of marginalisation in data-space enables our framework to provide objective measures for shape assignment in model-spaces as well. We note that the map between model-space and data-space used in these shape-calculations should be explicitly acknowledged and defined in order to mitigate unintentional mis-estimations of the weight of particular models in data-space during the change of variables. Finally, it is worth mentioning that such model-spaces almost always exist for scientific analyses, even if they are only implicitly invoked within $M_{0}$.
The most complete implementations of these BMS-based shape calculations occur when using $M_{0}$s that specify every physically appropriate $\bm{x}_{k}$. However, this approach may not always be theoretically possible nor computationally feasible if an effectively infinite number of $\bm{x}_{k}$s exist. In such situations, it is worth noting that, depending on the precision required by a particular analysis method, the full set of templates may not be required to obtain effective results. Importantly, a major benefit of BMS is that we can determine which $\bm{x}_{k}$ among a set of approximate templates best describes the shape of $\bm{y}$, even if none are ‘exact’. Additionally, the BMS expression in Eqn. (7) can be rearranged into a function of the log difference of evidence expressions (i.e., a Bayes’ factor) between a test $\bm{x}_{k}$ and an appropriate control $\bm{x}_{k}$ (or a ‘null’ model), which yields an effective cost function for the direct optimisation of a single $\bm{x}_{k}$ (see Supplemental Materials, Section 3). Overall, the most powerful aspect of the BMS-based shape calculations described here is that by considering different $M_{0}$s, an analysis method can be optimised for completeness (where all appropriate templates are enumerated) or efficiency (where only a test and a control template are considered), or for a trade-off between the two (using only a restricted set of templates), as the situation demands. This flexibility is a large reason why our framework can be effectively adapted to mimic nearly any of the subjective, expert-based analysis methods that it is meant to replace. Furthermore, the ability to easily disseminate the $\{\bm{x}\}$ used in an analysis means that methods can be readily shared, critiqued, and reproduced. Together, these especially powerful aspects of our framework make it extremely straightforward to implement tailor-made, shape-based analysis methods for new experiments and applications.