Declarative Modeling and Bayesian Inference of Dark Matter Halos

by   Gabriel Kronberger, et al.
FH Oberösterreich

Probabilistic programming allows specification of probabilistic models in a declarative manner. Recently, several new software systems and languages for probabilistic programming have been developed on the basis of newly developed and improved methods for approximate inference in probabilistic models. In this contribution a probabilistic model for an idealized dark matter localization problem is described. We first derive the probabilistic model for the inference of dark matter locations and masses, and then show how this model can be implemented using BUGS and Infer.NET, two software systems for probabilistic programming. Finally, the different capabilities of both systems are discussed. The presented dark matter model includes mainly non-conjugate factors, thus, it is difficult to implement this model with Infer.NET.


page 1

page 2

page 3

page 4


Deployable probabilistic programming

We propose design guidelines for a probabilistic programming facility su...

Inference Compilation and Universal Probabilistic Programming

We introduce a method for using deep neural networks to amortize the cos...

Bayesian additive regression trees for probabilistic programming

Bayesian additive regression trees (BART) is a non-parametric method to ...

When do Numbers Really Matter?

Common wisdom has it that small distinctions in the probabilities quanti...

Sound Probabilistic Inference via Guide Types

Probabilistic programming languages aim to describe and automate Bayesia...

Absolutely No Free Lunches!

This paper is concerned with learners who aim to learn patterns in infin...

A Probabilistic Programming Idiom for Active Knowledge Search

In this paper, we derive and implement a probabilistic programming idiom...

1 Introduction

Recently, there has been a growing interest in declarative probabilistic modeling which has led to the development of several software systems for probabilistic modeling and Bayesian inference such as Stan [13], FACTORIE [8], Infer.NET [9], or PRISM [12]

. Many of these systems provide a declarative modeling language for the definition of probabilistic models, which allows to define random variables and their relations in a way similar to computer programs. Thus, the term

probabilistic programming is frequently used to refer to the implementation of probabilistic models in such systems. The common idea is to implement the probabilistic model declaratively, without specifying how inference should be performed in the model. Instead the underlying inference engine is responsible for the execution of an appropriate inference algorithm and can potentially adapt the inference procedure to specific models (e.g., to improve accuracy of efficiency).

This approach to probabilistic modeling is not a recent idea; BUGS [5], a software system for Bayesian modeling and inference using Gibbs sampling, is already more than twenty years old [6]

and has become the de facto standard for probabilistic programming. BUGS defines its own modeling language, and relies on the fact, that Gibbs sampling is a very general Markov-chain Monte Carlo (MCMC) method and allows inference in a large class of probabilistic models. Thus, BUGS imposes almost no constraints on models and supports a large set of models, including analytically intractable models with non-conjugate or improper priors. However, MCMC methods often suffer from slow convergence especially for high-dimensional models or strongly correlated parameters. This drawback of Gibbs sampling has been a limiting factor for probabilistic programming with BUGS.

However, recent research results have led to several new and improved methods for approximate inference in probabilistic models, such as expectation propagation (EP) [10], variational message passing (VMP) [15], and improved sampling techniques including Hamiltonian Monte Carlo [11] and NUTS [2]. Several software systems have been developed, which incorporate these improved methods and can be used instead of BUGS.

In this contribution we discuss a Bayesian model for an idealized formulation of the dark matter localization problem, and show how this model can be implemented using BUGS and Infer.NET. The aim is to highlight and discuss the differences between BUGS and Infer.NET on the basis of a moderately complex model. A summary of the different capabilities of BUGS and Infer.NET, as well as a comparison of inference results, have also been given in [14].

1.1 Dark Matter Localization

A large fraction of the total mass in the universe is made up of so-called dark matter. Dark matter does not emit or absorb light but can be detected indirectly through its gravitational field. The existence and substance of dark matter is one of the unanswered questions of astrophysics, and a lot of effort is spent on improving methods to detect dark matter, and on studying its distribution in the universe. For instance in a recent publication a map of the distribution of dark mapper in the universe is discussed [7].

Dark matter can be detected through the gravitational lensing effect [4], which occurs because the gravitational force of large masses has a bending effect on light. Because of the gravitational lensing effect, objects behind a mass appear displaced to an observer and even multiple images of the same object might appear. Additionally, the apparent shape of larger objects such as galaxies is altered by the gravitational lensing effect.

The main aim of this contribution is to show, how a moderately complex probabilistic model, such as the dark matter localization model, can be implemented using software systems for probabilistic programming, in order to highlight and discuss the capabilities of such systems. We do not aim to derive a model that can be actually used for dark matter localization. However, it should be noted that e.g., LENSTOOL, a software system which has actually been used for calculating mass distribution profiles based on real images, also implements a Bayesian model and MCMC sampling [3].

1.2 Synthetic Data for Dark Matter Localization

We use a synthetic data set for the experiments presented in this contribution as we do not aim to improve on established models for dark matter localization. The data set has been generated for the “Observing Dark Worlds” competition hosted on Kaggle111Competition organizers: David Harvey and Thomas Kitching, Observing Dark Worlds Competition, by simulating dark matter halos, galaxies and the gravitational lensing effect. Distortions that would occur in real images e.g., through atmospheric effects or telescopic lenses, are ignored. For the purpose of the competition, real image data could not be used as it is necessary to compare solutions for dark matter halo locations. The data set is composed of 300 simulated skies and either one, two, or three dark matter halos. The data for each sky contains locations and ellipticities of between 300 and 740 galaxies. The ellipticity is specified using two components: the ellipticity along the x-axis , and the ellipticity along a 45-degree angle to the x-axis .

2 Model Formulation

In the following we describe the probabilistic graphical model for Bayesian inference of dark matter halo locations and masses from the observed locations and ellipticities of galaxies as specified in the synthetic data set.

It should be noted, that the model described below is the model that has been used by the author for the dark worlds competition. This model is very similar to the model used by the competition winner222Tim Salimans’ description of his model can be found at, but differs in relevant details. In particular, the model below uses different priors.

The goal in the dark matter localization problem is to determine probability distributions for halo locations

given observed galaxy locations and ellipticities

. Using Bayes’ theorem we can derive the posterior distribution of halo locations from the likelihood of observed locations and ellipticities times the prior. The mass of a halo determines the strength of the gravitational lensing effect. Therefore, the halo mass has to be included in the model as latent variable. This leads to the following model

There is no prior information about the locations of halos, so a flat uniform prior is used. For the halo mass a broad gamma prior is used: .

We assume that the galaxy locations are independent from halo locations and masses so the likelihood can be transformed:

The ellipticity is given as a vector with two components

which are independent, and the empirical distribution of ellipticities is very close to a zero mean normal distribution with variance

. Therefore, the likelihood function for observed ellipticities can be expressed as a normal likelihood with the mean given by the strength of the lensing effect and a constant variance .

We assume that the strength of the lensing effect can be modeled using a simple function, which increases linearly with the mass of the halo and inversely with the distance to the halo center.

The final model is shown as a probabilistic graphical model in Figure 1, using plate notation to represent galaxies and halos. In the graphical representation several details that are necessary for the implementation, such as the distance of galaxies from halos and the angle of the force vector, are not shown.

Figure 1: Probabilistic graphical model for the gravitational lensing effect that can be used to infer dark matter locations and masses. The strength of the lensing effect is a result of the masses and locations of halos and effects the observed ellipticities of galaxies.

For the implementation of the model the total force has to be allocated to the two components of ellipticity, which describe the elongation of the galaxy in the x-direction () and the elongation along a angle. The gravitational lensing effect leads to tangential elongation of the ellipticity. Thus, we map the tangential force to the two components by multiplying with and , respectively. is the angle of the vector from the halo to the galaxy.

2.1 Implementation using BUGS

It is rather straightforward to transform the probabilistic model to BUGS syntax. The model can be defined as follows:

  for( i in 1 : G ) {
    for( h in 1 : H ) {
      dx[i , h]    <- gx[i] - loc[h , 1]
      dy[i , h]    <- gy[i] - loc[h , 2]
      dist[i , h]  <- sqrt(dx[i , h] * dx[i , h] +
                           dy[i , h] * dy[i , h])
      phi[i , h]   <- atan2(dy[i , h], dx[i , h])
      iDist[i , h] <- 1.0 / dist[i , h]
      f[i , h]     <- mass[h] * iDist[i , h]
      f1[i , h]    <- -f[i , h] * cos(2 * phi[i , h])
      f2[i , h]    <- -f[i , h] * sin(2 * phi[i , h])
    mu1[i] <- sum(f1[i , ])
    mu2[i] <- sum(f2[i , ])
    e1[i] ~ dnorm(mu1[i], 0.05)
    e2[i] ~ dnorm(mu2[i], 0.05)
  for( h in 1 : H ) {
    mass[h]    ~ dgamma(0.001, 0.001)
    loc[h , 1] ~ dunif(0, 4200)
    loc[h , 2] ~ dunif(0, 4200)

The only additional steps that are necessary to infer and are loading initial values for all unobserved variables and loading the data for galaxy ellipticities and locations . However, convergence is very slow when sampling this model within BUGS. Additionally, BUGS does not support the atan2() function in models. Fortunately, the source code of BUGS is available, so it is possible to add support for this function rather easily.

2.2 Implementation using Infer.NET

Infer.NET333Infer.NET version 2.5 is available from
 [9] allows declarative specification of probabilistic models and provides EP [10] and VMP [15] for approximate inference. Probabilistic models can be implemented directly in C#. The model is transformed transparently by the Infer.NET compiler to C# source code, which is then compiled to CLR byte code using the C# compiler. Compared to BUGS, it is much easier to use such models from existing code, as long as the application is based on the .NET platform. Additionally, inference is fast because the code for model inference is compiled. The model can be implemented in the following way:

// not shown: variable declarations
// [...]
sigma = Variable.New<double>();
evidence = Variable.Bernoulli(0.5);

IfBlock block = Variable.If(evidence);
using (Variable.ForEach(g)) {
  using (Variable.ForEach(h)) {
    // factors for the following are not implemented
    // dx[g][h].SetTo(g_x[g] - loc_x[h]);
    // dy[g][h].SetTo(g_y[g] - loc_y[h]);
    // phi.SetTo(Variable.Atan2(dy[g][h], dx[g][h]));
    // cos2phi[g][h].SetTo(Variable.Cos(2 * phi[g][h]));
    // sin2phi[g][h].SetTo(Variable.Sin(2 * phi[g][h]));
    // invDist[g][h].SetTo(1.0 / Variable.Sqrt(dx[g][h] * dx[g][h] +
    //                                         dy[g][h] * dy[g][h]));
    f1[g][h].SetTo(-(mass[h] * invDist[g][h] * cos2phi[g][h]));
    f2[g][h].SetTo(-(mass[h] * invDist[g][h] * sin2phi[g][h]));
   Variable.GaussianFromMeanAndPrecision(Variable.Sum(f1[g]), 20));
   Variable.GaussianFromMeanAndPrecision(Variable.Sum(f2[g]), 20));

Similarly to the BUGS implementation, arrays of random variables are used for galaxy locations and ellipticities, as well as for halo locations and masses. Looping over ranges can be accomplished with the Variable.ForEach factor; in the example two variants to handle blocks in Infer.NET are shown. Either the block is manually opened and closed, as shown for the IfBlock, or the using-syntax is used to manage blocks.

EP and VMP are able to exploit regularities in the model, so that inference can be performed much faster than would be possible with MCMC techniques. EP and VMP allow inference also for large scale models, such as document topic modeling with latent Dirichlet allocation. The drawback of both methods is that they are much less general than e.g., Gibbs sampling. Thus, the set of probabilistic models that can be used in Infer.NET is rather constrained [14]. However, in contrast to BUGS, Infer.NET also supports conditional blocks.

Even tough EP can potentially also be used to infer marginals when the factors are non-conjugate, this is in general not very efficient. So, Infer.NET typically does not contain such factors. Implementing custom factors for EP inference is rather difficult.

The model uses mainly non-conjugate factors. Thus, it is not possible to perform inference using the original model with the standard installation of Infer.NET. It would be necessary to implement new distribution types, factors and message operators to support inference for this model. Instead, the model is simplified, so that it basically represents a simple likelihood function. Infer.NET is only used to infer the likelihood (evidence). The parameters of the model, and are optimized w.r.t. likelihood using the CMA-ES optimization algorithm [1]. Figure 2 shows locations and shapes of galaxies and the actual location of the dark matter halos for two skies. The predicted halo locations are marked with a green circle.

Figure 2: Two simulated skies showing galaxies and the actual (red cross) and predicted (green circle) center location of the dark matter halos.

3 Summary and Discussion

In this contribution a Bayesian model for gravitational lensing is described that the author used for the “Observing Dark Worlds” competition. The model is similar to the winning model but uses e.g., different priors. The gravitational lensing model can be easily implemented in BUGS. Only the function is missing, however, it is easy to add support for this function to the open-source version of BUGS. Sampling converges only slowly for this model in BUGS.

Infer.NET uses EP or VMP for approximate inference, and the models can be implemented directly e.g., in C#. Inference with EP and VMP is very efficient for certain types of models, but in general less flexible than e.g., Gibbs sampling. Because of these restrictions, it would be necessary to implement custom factors and distributions, to implement the full gravitational lensing model using Infer.NET. Instead, we performed simple ML optimization of dark halo locations and masses, based on a simplified implementation of the model in Infer.NET, using CMA-ES as optimization algorithm. The results are competitive, taking a spot in the top 10% of submitted solutions.

In this contribution we have only discussed the implementation of the model using BUGS and Infer.NET. However, several other software systems for probabilistic programming, e.g. Stan or FACTORIE, can also be used instead. It would certainly be interesting to implement the model also in these systems, and to compare the results and performance. Additionally, it would be interesting to discuss the implementation of factors, which are necessary for inference with EP, in more detail.


  • [1]

    Hansen, N., Ostermeier, A.: Completely derandomized self-adaptation in evolution strategies. Evolutionary computation 9(2), 159–195 (2001)

  • [2] Hoffman, M.D., Gelman, A.: The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo. ArXiv e-prints (Nov 2011)
  • [3] Jullo, E., Kneib, J.P., Limousin, M., Eliasdottir, A., Marshall, P., Verdugo, T.: A Bayesian approach to strong lensing modelling of galaxy clusters. New Journal of Physics 9(447) (2007)
  • [4] Kaiser, N., Squires, G.: Mapping the dark matter with weak gravitational lensing. Astrophysical Journal 404(2), 441–450 (1993)
  • [5] Lunn, D., Jackson, C., Best, N., Spiegelhalter, D.J., Thomas, A.: The BUGS book: A practical introduction to Bayesian analysis, vol. 98. Chapman & Hall (2012)
  • [6] Lunn, D., Spiegelhalter, D., Thomas, A., Best, N.: The BUGS project: Evolution, critique and future directions. Statistics in medicine 28(25), 3049–3067 (2009)
  • [7] Massey, R., Rhodes, J., Ellis, R., Scoville, N., Leauthaud, A., Finoguenov, A., Capak, P., Bacon, D., Aussel, H., Kneib, J.P., Koekemoer, A., McCracken, H., Mobasher, B., Pires, S., Refregier, A., Sasaki, S., Starck, J.L., Taniguchi, Y., Taylor, A., Taylor, J.: Dark matter maps reveal cosmic scaffolding. Nature 445, 286–290 (2007)
  • [8] McCallum, A., Schultz, K., Singh, S.: FACTORIE: Probabilistic programming via imperatively defined factor graphs. In: Neural Information Processing Systems (NIPS) (2009)
  • [9] Minka, T., Winn, J., Guiver, J., Knowles, D.: Infer.NET 2.5 (2012), Microsoft Research Cambridge.
  • [10]

    Minka, T.P.: Expectation propagation for approximate Bayesian inference. In: Proceedings of the Seventeenth conference on Uncertainty in artificial intelligence. pp. 362–369. Morgan Kaufmann Publishers Inc. (2001)

  • [11] Neal, R.M.: MCMC using Hamiltonian dynamics. In: Brooks, S., Gelman, A., Jones, G.L., Meng, X.L. (eds.) Handbook of Markov Chain Monte Carlo, pp. 113–162. Chapman and Hall/CRC (2011)
  • [12]

    Sato, T., Kameya, Y.: New advances in logic-based probabilistic modeling by PRISM. In: Raedt, L., Frasconi, P., Kersting, K., Muggleton, S. (eds.) Probabilistic Inductive Logic Programming, Lecture Notes in Computer Science, vol. 4911, pp. 118–155. Springer Berlin Heidelberg (2008)

  • [13] Stan Development Team: Stan: A C++ Library for Probability and Sampling, Version 1.3 (2013),
  • [14] Wang, S., Wand, M.: Using infer.NET for statistical analyses. The American Statistician 65(2), 115–126 (2011)
  • [15]

    Winn, J., Bishop, C.M.: Variational message passing. Journal of Machine Learning Research 6(1), 661 (2006)