Noise models

In socca, the likelihood function to be used during the inference process is fully specified by the noise model associated with the input data. This encodes the assumed statistical properties of the noise affecting the observations and therefore plays a central role in defining how model predictions are compared to data themselves. As shown in “Getting started”, any noise model should be linked to the input data by passing the noise object to the corresponding socca.data.Image instance:

>>> img = Image(..., noise=noise)

At present, socca provides a limited but representative set of Gaussian noise models, covering the most common use cases encountered in astronomical imaging and map-based analyses. Below is a brief description of the noise models currently implemented in socca.

Note

Several noise models are under active development and testing, and will be added to socca in future releases. These include, among others, Poisson noise models for photon-counting data (e.g., X-ray observations).

Important

If the noise model is not explicitly specified when instantiating a socca.data.Image object, socca will use the Normal noise model by default and compute the per-pixel standard deviation from the median absolute deviation of all pixels in the input data that are not masked out.

Normal: uncorrelated normal noise

The Normal noise model implements a multivariate normal likelihood, assuming that each spatial pixel in the input data is statistically independent from all others and affected by independent Gaussian noise with zero mean and known standard deviation. This is the simplest noise model available in socca and is appropriate when pixel-to-pixel correlations can be neglected and the noise variance is either uniform or known on a per-pixel basis.

The Normal noise model can be instantiated without any arguments. In such a case, the per-pixel standard deviation is automatically estimated from the median absolute deviation of all pixels in the input data, under the assumption that the noise is homogeneous across the image.

>>> from socca.noise import Normal
>>> noise = Normal()

Alternatively, the noise amplitude can be specified explicitly using any of the following equivalent keyword arguments. Internally, all parametrizations are converted to a standard-deviation map.

Argument	Aliases	Meaning	Relation
sigma	sig, std, rms, stddev	Standard deviation	σ
variance	var	Variance	σ²
weight	wht, wgt, weights, invvar	Inverse variance	1 / σ²

These can be provided either as scalar values (assuming homogeneous noise levels) or as arrays with the same shape as the input data (allowing for spatially varying noise properties). If a string is provided for any of these arguments, it is interpreted as the path to a FITS file containing the corresponding maps.

NormalCorrelated: correlated normal noise

The NormalCorrelated noise model generalizes the uncorrelated Gaussian case by allowing for for spatial correlations between different pixels. The noise correlations are described explicitly through a full covariance matrix defined on the flattened data vector.

In this case, information about the covariance structure of the noise must always be provided explicitly when instantiating the model. This can be done by providing the full covariance matrix via the cov argument, or directly its inverse via the icov argument. In the former case, the inverse covariance matrix is computed internally using jax.numpy.linalg.inv. The inverse covariance matrix is used directly when evaluating the (normalized) likelihood.

>>> from socca.noise import NormalCorrelated
>>>
>>> # From a covariance matrix
>>> noise_from_cov = NormalCorrelated(cov=covariance_matrix)
>>>
>>> # From an inverse covariance matrix
>>> noise_from_inv = NormalCorrelated(icov=inverse_covariance_matrix)

Both the cov and icov arguments should be provided as 2D jax.numpy arrays of shape (nx*ny,nx*ny), where (nx,ny) is the shape of the input data. The two options are equivalent, and the choice is left to the user depending on which representation is more convenient to obtain.

Alternatively, NormalCorrelated offers the possibility of estimating the noise covariance matrix internally from a set of noise realizations provided to socca via the cube argument.

>>> from socca.noise import NormalCorrelated
>>> noise = NormalCorrelated(cube=noise_realizations)

Here, noise_realizations should be a 3D jax.numpy array of shape (nr,nx,ny), where nr is the number of independent noise maps available. Once computed, the covariance matrix is smoothed using a simple convolutional kernel to reduce the impact of noise fluctuations. The amount of smoothing can be controlled via the smooth argument, denoting the number of iterations of the smoothing operation. By default, it is performed 3 times, using a simple 5-point stencil kernel. A custom smoothing kernel can also be provided via the kernel argument.

Warning

This noise model explicitly evaluates the likelihood using the full inverse covariance matrix and, therefore, is generally computationally more expensive than the uncorrelated Normal noise model.

NormalFourier: uncorrelated Fourier-space normal noise

The NormalFourier noise model implements a Gaussian likelihood assuming that the noise is independent in Fourier space. Under this assumption, the noise covariance is diagonal in Fourier space and is defined directly on the Fourier coefficients of the data. This noise model is appropriate for approximately stationary noise processes, for which correlations in real space become simple in Fourier space. Typical use cases include map-based data originating from time-ordered or scan-based observations.

As in the case of NormalCorrelated, the noise properties can be specified either by providing the noise covariance in Fourier space via the cov argument, by providing its inverse via the icov argument, or by providing a set of noise realizations via the cube argument, from which the Fourier-space covariance is estimated internally. When a noise realization cube is provided, each realization is Fourier-transformed (optionally after apodization), and the covariance is computed as the mean squared modulus of the Fourier coefficients. The resulting covariance can optionally be smoothed before inversion.

>>> from socca.noise import NormalFourier
>>>
>>> # From a Fourier-space covariance
>>> noise = NormalFourier(cov=fourier_covariance)
>>>
>>> # From an inverse Fourier-space covariance
>>> noise = NormalFourier(icov=inverse_fourier_covariance)
>>>
>>> # From a cube of noise realizations
>>> noise = NormalFourier(cube=noise_cube, ftype="real")

The optional argument ftype controls the type of Fourier transform used when estimating the covariance from the noise cube. By default, a real-to-complex 2D Fourier transform is used (ftype="real"), in order to minimize the memory usage and computational cost. If needed, it is possible to use a complex-to-complex 2D Fourier transform by specifying ftype="full".

Other optional keyword arguments include apod, which specifies an apodization map applied to the data before Fourier transforming (though, by default, no apodization is applied). Finally, as in the case of NormalCorrelated, the smooth and kernel arguments control the optional smoothing of the Fourier-space covariance matrix when estimated from the noise realizations. As above, the default is 3 smoothing iterations using a 5-point stencil kernel.

Warning

The likelihood is evaluated by Fourier transforming the apodized residuals between model and data, weighting each Fourier mode by the inverse noise covariance, and summing over the Fourier modes with non-zero inverse covariance. Given the Fourier-space nature of this noise model, it is not possible to handle masked pixels in the input data.

NormalRI: radio-interferometric noise

The NormalRI noise model implements an image-space approximation to the Fourier-space likelihood for radio-interferometric data. Rather than working directly in the visibility domain, this model uses the dirty beam (i.e., the interferometric response function) to approximate the noise covariance structure in image space. This approach enables efficient likelihood evaluation while accounting for the correlated noise introduced by the interferometric imaging process. The specific implementation is based on the derivations by Powell et al. (2020) and Zhang et al. (2025), to which we refer to for an extensive discussion of the theoretical framework.

Important

When using the NormalRI noise model, the input data should be a naturally-weighted dirty image (i.e., without any deconvolution applied) and should not be corrected for the primary beam (i.e., the antenna pattern) response. If the primary beam is expected to have a significant effect on the signal, a model should be provided as the response attribute of the input socca.data.Image object.

The noise amplitude can be specified using the same keyword arguments as the Normal noise model (see above), but does not support image-like inputs. If no noise amplitude is provided, the per-pixel standard deviation is estimated from the median absolute deviation of the unmasked pixels.

>>> from socca.noise import NormalRI
>>>
>>> # With explicit noise level
>>> noise = NormalRI(sigma=0.1)
>>>
>>> # With automatic noise estimation
>>> noise = NormalRI()

Warning

The NormalRI noise model is designed for map-based analyses of radio-interferometric data, where the dirty beam encodes the sampling of the Fourier plane by the interferometer. It provides a computationally efficient alternative to full visibility-domain fitting, at the cost of an approximate treatment of the noise correlations. Current tests indicate that this approach is accurate enough for sources that are not too extended compared to the maximum recoverable scale in the specific interferometric dataset. For more extended sources, the approximation may break down and a full visibility-domain analysis may be necessary.

Attention

Some tests have exposed some stability issues in the case of priors that allow for sources with characteristic scales much larger than the maximum recoverable scale and the field of view of the interferometric dataset. In such cases, the sampler might explore regions of the parameter space where the model predictions are dominated by large-scale features that are not well constrained by the data, leading to an unstable behavior of the sampling process. This isssue is currently under investigation and will be addressed in future releases of socca. In the meantime, we recommend using more informative priors that restrict the parameter space to regions where the model predictions are reasonably constrained by the data, especially when dealing with sources that are expected to be extended compared to the maximum recoverable scale.