Hamiltonian Monte Carlo Sampling Research Articles

On the effectiveness of partially deterministic Bayesian neural networks

AbstractBayesian neural networks (BNNs) with computationally expensive Hamiltonian Monte Carlo sampling methods are often considered to provide better predictive performance than the maximum a posterior (MAP) solution. Here, as an alternative to sampling all parameters of a BNN (full-random), we experimentally evaluate partially deterministic BNNs that fix some part of the neural network parameters to their MAP solution. In particular, we consider various strategies for fixing half, or all parameters of a layer to the MAP-solution. Over a wide variety of regression and classification tasks, we find that partially deterministic BNNs often significantly improve predictive performance over the MAP-solution, with up to around 24% reduction in negative log-likelihood. Notably, we also find that partially deterministic BNNs that fix half of the parameters in each layer can also reduce under-fitting of full-random BNNs, resulting in up to 7% reduction in negative log-likelihood.

A failure-dependence related stochastic crack growth modeling approach of competing cracking mode

Competing crack-induced fatigue fracture has been a typical failure mode especially in the very high cycle fatigue regime. However, the stochastic crack growth modeling related to failure-dependence has not been fully investigated. Here, a failure-dependence related stochastic crack growth modeling approach for competing cracking mode is proposed to address this issue. Firstly, copulas are used to model the dependent competing relationships among multiple fatigue fracture modes. The reliability analysis of multiple fatigue fracture modes is conducted from a copula perspective. Besides, the fatigue crack growth is modeled based on a nonlinear Wiener process, with unit-to-unit variability addressed by introducing random effects. Marginal reliability expressions are derived based on the Wiener process. Furthermore, a two-stage Bayesian inference method based on Hamiltonian Monte Carlo sampling is proposed to estimate model parameters. A Monte Carlo simulation study is conducted to validate the accuracy and robustness of the proposed inference method. Finally, the effectiveness of the proposed approach is proven through a real case study. It turns out that the ignorance of the competing relationships among multiple cracking modes leads to an underestimation of overall reliability. The accuracy of the model can be further improved with random effects considered.

Bayesian Modeling of Travel Times on the Example of Food Delivery: Part 2—Model Creation and Handling Uncertainty

The e-commerce sector is in a constant state of growth and evolution, particularly within its subdomain of online food delivery. As such, ensuring customer satisfaction is critical for companies working in this field. One way to achieve this is by providing an accurate delivery time estimation. While companies can track couriers via GPS, they often lack real-time data on traffic and road conditions, complicating delivery time predictions. To address this, a range of statistical and machine learning techniques are employed, including neural networks and specialized expert systems, with different degrees of success. One issue with neural networks and machine learning models is their heavy dependence on vast, high-quality data. To mitigate this issue, we propose two Bayesian generalized linear models to predict the time of delivery. Utilizing a linear combination of predictor variables, we generate a practical range of outputs with the Hamiltonian Monte Carlo sampling method. These models offer a balance of generality and adaptability, allowing for tuning with expert knowledge. They were compared with the PSIS-LOO criteria and WAIC. The results show that both models accurately estimated delivery times from the dataset while maintaining numerical stability. A model with more predictor variables proved to be more accurate.

Reliability analysis of complex systems using subset simulations with Hamiltonian Neural Networks

We present a new Subset Simulation approach using Hamiltonian neural network-based Monte Carlo sampling for reliability analysis. The proposed strategy combines the superior sampling of the Hamiltonian Monte Carlo method with computationally efficient gradient evaluations using Hamiltonian neural networks. This combination is especially advantageous because the neural network architecture conserves the Hamiltonian, which defines the acceptance criteria of the Hamiltonian Monte Carlo sampler. Hence, this strategy achieves high acceptance rates at low computational cost. Our approach estimates small failure probabilities using Subset Simulations. However, in low-probability sample regions, the gradient evaluation is particularly challenging. The remarkable accuracy of the proposed strategy is demonstrated on different reliability problems, and its efficiency is compared to the traditional Hamiltonian Monte Carlo method. We note that this approach can reach its limitations for gradient estimations in low-probability regions of complex and high-dimensional distributions. Thus, we propose techniques to improve gradient prediction in these particular situations and enable accurate estimations of the probability of failure. The highlight of this study is the reliability analysis of a system whose parameter distributions must be inferred with Bayesian inference problems. In such a case, the Hamiltonian Monte Carlo method requires a full model evaluation for each gradient evaluation and, therefore, comes at a very high cost. However, using Hamiltonian neural networks in this framework replaces the expensive model evaluation, resulting in tremendous improvements in computational efficiency.

Scalable gradients enable Hamiltonian Monte Carlo sampling for phylodynamic inference under episodic birth-death-sampling models.

Birth-death models play a key role in phylodynamic analysis for their interpretation in terms of key epidemiological parameters. In particular, models with piecewise-constant rates varying at different epochs in time, to which we refer as episodic birth-death-sampling (EBDS) models, are valuable for their reflection of changing transmission dynamics over time. A challenge, however, that persists with current time-varying model inference procedures is their lack of computational efficiency. This limitation hinders the full utilization of these models in large-scale phylodynamic analyses, especially when dealing with high-dimensional parameter vectors that exhibit strong correlations. We present here a linear-time algorithm to compute the gradient of the birth-death model sampling density with respect to all time-varying parameters, and we implement this algorithm within a gradient-based Hamiltonian Monte Carlo (HMC) sampler to alleviate the computational burden of conducting inference under a wide variety of structures of, as well as priors for, EBDS processes. We assess this approach using three different real world data examples, including the HIV epidemic in Odesa, Ukraine, seasonal influenza A/H3N2 virus dynamics in New York state, America, and Ebola outbreak in West Africa. HMC sampling exhibits a substantial efficiency boost, delivering a 10- to 200-fold increase in minimum effective sample size per unit-time, in comparison to a Metropolis-Hastings-based approach. Additionally, we show the robustness of our implementation in both allowing for flexible prior choices and in modeling the transmission dynamics of various pathogens by accurately capturing the changing trend of viral effective reproductive number.

Data-driven predictions of shield attitudes using Bayesian machine learning

Shield attitude prediction is essential in shield tunnelling to prevent excessive attitude deviations and is highly correlated with shield operational parameters. Furthermore, geological and geotechnical uncertainty are the primary sources of attitude deviation risks; however, these uncertainties were neglected in previous studies. In this paper, a novel data-driven Bayesian machine learning method for site-specific shield attitude predictions is proposed by simultaneously considering the spatial variability of geotechnical parameters and the generic cross-correlations between shield attitude parameters and operational parameters. The Gaussian process model is first used to effectively capture the spatial variability of geotechnical parameters using field borehole data. The multivariate probability distribution model is then constructed to characterise the generic cross-correlation between four attitude parameters, four geotechnical parameters and eight shield operational parameters. In the Bayesian theory framework, a hybrid strategy is proposed to fuse the generic data from other shield tunnelling projects with the ongoing site-specific data, ensuring that the shield attitude prediction is dominated by site-specific data when those data are abundant or by generic data when the site-specific data are sparse. Finally, the Hamiltonian Monte Carlo sampling method is employed to generate the posterior distributions of shield attitudes conditioned on geotechnical parameters and shield operational parameters. The proposed method has the benefit of rigorously quantifying the uncertainties of shield attitude predictions. The proposed approach is examined on a real-world tunnelling project in Nanjing City, indicating that the proposed method outperforms existing techniques in terms of both accuracy and robustness. Parametric analyses are conducted to explore the impacts of key factors on the performance of the proposed approach.

Detailed statistical analysis plan for a secondary Bayesian analysis of the SafeBoosC-III trial: a multinational, randomised clinical trial assessing treatment guided by cerebral oximetry monitoring versus usual care in extremely preterm infants

BackgroundExtremely preterm infants have a high mortality and morbidity. Here, we present a statistical analysis plan for secondary Bayesian analyses of the pragmatic, sufficiently powered multinational, trial—SafeBoosC III—evaluating the benefits and harms of cerebral oximetry monitoring plus a treatment guideline versus usual care for such infants.MethodsThe SafeBoosC-III trial is an investigator-initiated, open-label, randomised, multinational, pragmatic, phase III clinical trial with a parallel-group design. The trial randomised 1601 infants, and the frequentist analyses were published in April 2023. The primary outcome is a dichotomous composite outcome of death or severe brain injury. The exploratory outcomes are major neonatal morbidities associated with neurodevelopmental impairment later in life: (1) bronchopulmonary dysplasia; (2) retinopathy of prematurity; (3) late-onset sepsis; (4) necrotising enterocolitis; and (5) number of major neonatal morbidities (count of bronchopulmonary dysplasia, retinopathy of prematurity, and severe brain injury). The primary Bayesian analyses will use non-informed priors including all plausible effects. The models will use a Hamiltonian Monte Carlo sampler with 1 chain, a sampling of 10,000, and at least 25,000 iterations for the burn-in period. In Bayesian statistics, such analyses are referred to as ‘posteriors’ and will be presented as point estimates with 95% credibility intervals (CrIs), encompassing the most probable results based on the data, model, and priors selected. The results will be presented as probability of any benefit or any harm, Bayes factor, and the probability of clinical important benefit or harm. Two statisticians will analyse the blinded data independently following this protocol.DiscussionThis statistical analysis plan presents a secondary Bayesian analysis of the SafeBoosC-III trial. The analysis and the final manuscript will be carried out and written after we publicise the primary frequentist trial report. Thus, we can interpret the findings from both the frequentists and Bayesian perspective. This approach should provide a better foundation for interpreting of our findings.Trial registrationClinicalTrials.org, NCT03770741. Registered on 10 December 2018.

Fast reconstruction algorithm based on HMC sampling

In Noisy Intermediate-Scale Quantum (NISQ) era, the scarcity of qubit resources has prevented many quantum algorithms from being implemented on quantum devices. Circuit cutting technology has greatly alleviated this problem, which allows us to run larger quantum circuits on real quantum machines with currently limited qubit resources at the cost of additional classical overhead. However, the classical overhead of circuit cutting grows exponentially with the number of cuts and qubits, and the excessive postprocessing overhead makes it difficult to apply circuit cutting to large scale circuits. In this paper, we propose a fast reconstruction algorithm based on Hamiltonian Monte Carlo (HMC) sampling, which samples the high probability solutions by Hamiltonian dynamics from state space with dimension growing exponentially with qubit. Our algorithm avoids excessive computation when reconstructing the original circuit probability distribution, and greatly reduces the circuit cutting post-processing overhead. The improvement is crucial for expanding of circuit cutting to a larger scale on NISQ devices.

Hierarchical Bayesian pharmacometrics analysis of Baclofen for alcohol use disorder

Alcohol use disorder (AUD), also called alcohol dependence, is a major public health problem, affecting almost of the world’s population. Baclofen, as a selective receptor agonist, has emerged as a promising drug for the treatment of AUD. However, the inter-trial, inter-individual and residual variability in drug concentration over time in a population of patients with AUD is unknown. In this study, we use a hierarchical Bayesian workflow to estimate the parameters of a pharmacokinetic (PK) population model from Baclofen administration to patients with AUD. By monitoring various convergence diagnostics, the probabilistic methodology is first validated on synthetic longitudinal datasets and then applied to infer the PK model parameters based on the clinical data that were retrospectively collected from outpatients treated with oral Baclofen. We show that state-of-the-art advances in automatic Bayesian inference using self-tuning Hamiltonian Monte Carlo (HMC) algorithms provide accurate and decisive predictions on Baclofen plasma concentration at both individual and group levels. Importantly, leveraging the information in prior provides faster computation, better convergence diagnostics, and substantially higher out-of-sample prediction accuracy. Moreover, the root mean squared error as a measure of within-sample predictive accuracy can be misleading for model evaluation, whereas the fully Bayesian information criteria correctly select the true data generating parameters. This study points out the capability of non-parametric Bayesian estimation using adaptive HMC sampling methods for easy and reliable estimation in clinical settings to optimize dosing regimens and efficiently treat AUD.

Scalable Bayesian Divergence Time Estimation With Ratio Transformations.

Divergence time estimation is crucial to provide temporal signals for dating biologically important events from species divergence to viral transmissions in space and time. With the advent of high-throughput sequencing, recent Bayesian phylogenetic studies have analyzed hundreds to thousands of sequences. Such large-scale analyses challenge divergence time reconstruction by requiring inference on highly correlated internal node heights that often become computationally infeasible. To overcome this limitation, we explore a ratio transformation that maps the original $N-1$ internal node heights into a space of one height parameter and $N-2$ ratio parameters. To make the analyses scalable, we develop a collection of linear-time algorithms to compute the gradient and Jacobian-associated terms of the log-likelihood with respect to these ratios. We then apply Hamiltonian Monte Carlo sampling with the ratio transform in a Bayesian framework to learn the divergence times in 4 pathogenic viruses (West Nile virus, rabies virus, Lassa virus, and Ebola virus) and the coralline red algae. Our method both resolves a mixing issue in the West Nile virus example and improves inference efficiency by at least 5-fold for the Lassa and rabies virus examples as well as for the algae example. Our method now also makes it computationally feasible to incorporate mixed-effects molecular clock models for the Ebola virus example, confirms the findings from the original study, and reveals clearer multimodal distributions of the divergence times of some clades of interest.

Consistency tests of field level inference with the EFT likelihood

Analyzing the clustering of galaxies at the field level in principle promises access to all the cosmological information available. Given this incentive, in this paper we investigate the performance of field-based forward modeling approach to galaxy clustering using the effective field theory (EFT) framework of large-scale structure (LSS). We do so by applying this formalism to a set of consistency and convergence tests on synthetic datasets. We explore the high-dimensional joint posterior of LSS initial conditions by combining Hamiltonian Monte Carlo sampling for the field of initial conditions, and slice sampling for cosmology and model parameters. We adopt the Lagrangian perturbation theory forward model from [1], up to second order, for the forward model of biased tracers. We specifically include model mis-specifications in our synthetic datasets within the EFT framework.We achieve this by generating synthetic data at a higher cutoff scale Λ0, which controls which Fourier modes enter the EFT likelihood evaluation, than the cutoff Λ used in the inference. In the presence of model mis-specifications, we find that the EFT framework still allows for robust, unbiased joint inference of a) cosmological parameters — specifically, the scaling amplitude of the initial conditions — b) the initial conditions themselves, and c) the bias and noise parameters. In addition, we show that in the purely linear case, where the posterior is analytically tractable, our samplers fully explore the posterior surface. We also demonstrate convergence in the cases of nonlinear forward models.Our findings serve as a confirmation of the EFT field-based forward model framework developed in [2-7], and as another step towards field-level cosmological analyses of real galaxy surveys.

Learning, inference, and prediction on probability density functions with constrained Gaussian processes

Probability density functions (PDFs) play an important role in machine learning, statistics, and in many real-world applications. In this paper, we introduce a new framework to learn, infer and predict nonparametric PDFs with a Gaussian process prior. This is a challenging problem since there is no explicit formulas for conditional expectations and covariances. The proposed methods have two main advantages: Analyzing PDFs with an intrinsic manifold structure and establishing the connection with the Hilbert sphere when endowed with an appropriate metric. Therefore, all solutions are valid PDFs with non-negative values and integrate to one. We formulate the problem as constrained spherical Gaussian processes leading to an efficient solution with a spherical Hamiltonian Monte Carlo (HMC) sampling. We test and evaluate different strategies with extensive experiments on both simulations and real dataset.

Analysis of the Hamiltonian Monte Carlo genotyping algorithm on PROVEDIt mixtures including a novel precision benchmark.

We provide an internal validation study of a recently published precise DNA mixture algorithm based on Hamiltonian Monte Carlo sampling (Susik etal., 2022). We provide results for all 428 mixtures analysed by Riman etal. (2021) and compare the results with two state-of-the-art software products: STRmix™v2.6 and Euroformix v3.4.0. The comparison shows that the Hamiltonian Monte Carlo method provides reliable values of likelihood ratios (LRs) close to the other methods. We further propose a novel large-scale precision benchmark and quantify the precision of the Hamiltonian Monte Carlo method, indicating its improvements over existing solutions. Finally, we analyse the influence of the factors discussed by Buckleton etal. (2022).

Modeling Contingent Decision Behavior: A Bayesian Nonparametric Preference-Learning Approach

We propose a preference-learning algorithm for uncovering Decision Makers’ (DMs’) contingent evaluation strategies in the context of multiple criteria sorting. We assume the preference information in the form of holistic assignment examples derived from the analysis of alternatives’ performance vectors and textual descriptions. We characterize the decision policies using a mixture of threshold-based, value-driven preference models and associated latent topics. The latter serve as the stimuli underlying the contingency in decision behavior. Such a probabilistic model is constructed by using a flexible and nonparametric Bayesian framework. The proposed method adopts a hierarchical Dirichlet process as the prior so that a group of DMs can share a countably infinite number of contingent models and topics. For all DMs, it automatically identifies the components representing their evaluation strategies adequately. The posterior is summarized by using the Hamiltonian Monte Carlo sampling method. We demonstrate the method’s practical usefulness in a real-world recruitment problem considered by a Chinese IT company. We also compare the approach with counterparts that use a single preference model, implement the parametric framework, or consider each DM’s preferences individually. The results indicate that our approach performs favorably in both interpreting DMs’ contingent decision behavior and recommending decisions on new alternatives. Furthermore, the approach’s performance and robustness are investigated through a computational experiment involving real-world data sets. History: Accepted by Ram Ramesh, Area Editor for Data Science and Machine Learning. Funding: J. Liu received financial support from the National Natural Science Foundation of China [Grants 72071155 and 71701160]. M. Kadziński received financial support from the Polish National Science Center under the SONATA BIS project [Grant DEC-2019/34/E/HS4/00045]. X. Liao received financial support from the National Natural Science Foundation of China [Grant 71872144]. Supplemental Material: The software that supports the findings of this study is available within the paper and its Supplemental Information ( https://pubsonline.informs.org/doi/suppl/10.1287/ijoc.2023.1292 ) as well as from the IJOC GitHub software repository ( https://github.com/INFORMSJoC/2021.0328 ) at ( http://dx.doi.org/10.5281/zenodo.7608750 ).

Variational Hybrid Monte Carlo for Efficient Multi-Modal Data Sampling

The Hamiltonian Monte Carlo (HMC) sampling algorithm exploits Hamiltonian dynamics to construct efficient Markov Chain Monte Carlo (MCMC), which has become increasingly popular in machine learning and statistics. Since HMC uses the gradient information of the target distribution, it can explore the state space much more efficiently than random-walk proposals, but may suffer from high autocorrelation. In this paper, we propose Langevin Hamiltonian Monte Carlo (LHMC) to reduce the autocorrelation of the samples. Probabilistic inference involving multi-modal distributions is very difficult for dynamics-based MCMC samplers, which is easily trapped in the mode far away from other modes. To tackle this issue, we further propose a variational hybrid Monte Carlo (VHMC) which uses a variational distribution to explore the phase space and find new modes, and it is capable of sampling from multi-modal distributions effectively. A formal proof is provided that shows that the proposed method can converge to target distributions. Both synthetic and real datasets are used to evaluate its properties and performance. The experimental results verify the theory and show superior performance in multi-modal sampling.

A monotone single index model for missing-at-random longitudinal proportion data

Beta distributions are commonly used to model proportion valued response variables, often encountered in longitudinal studies. In this article, we develop semi-parametric Beta regression models for proportion valued responses, where the aggregate covariate effect is summarized and flexibly modeled, using a interpretable monotone time-varying single index transform of a linear combination of the potential covariates. We utilize the potential of single index models, which are effective dimension reduction tools and accommodate link function misspecification in generalized linear mixed models. Our Bayesian methodology incorporates the missing-at-random feature of the proportion response and utilize Hamiltonian Monte Carlo sampling to conduct inference. We explore finite-sample frequentist properties of our estimates and assess the robustness via detailed simulation studies. Finally, we illustrate our methodology via application to a motivating longitudinal dataset on obesity research recording proportion body fat.

Almanac: Weak Lensing power spectra and map inference on the masked sphere

We present a field-based signal extraction of weak lensing from noisy observations on the curved and masked sky. We test the analysis on a simulated Euclid-like survey, using a Euclid-like mask and noise level. To make optimal use of the information available in such a galaxy survey, we present a Bayesian method for inferring the angular power spectra of the weak lensing fields, together with an inference of the noise-cleaned tomographic weak lensing shear and convergence (projected mass) maps. The latter can be used for field-level inference with the aim of extracting cosmological parameter information including non-gaussianity of cosmic fields. We jointly infer all-sky $E$-mode and $B$-mode tomographic auto- and cross-power spectra from the masked sky, and potentially parity-violating $EB$-mode power spectra, up to a maximum multipole of $\ell_{\rm max}=2048$. We use Hamiltonian Monte Carlo sampling, inferring simultaneously the power spectra and denoised maps with a total of $\sim 16.8$ million free parameters. The main output and natural outcome is the set of samples of the posterior, which does not suffer from leakage of power from $E$ to $B$ unless reduced to point estimates. However, such point estimates of the power spectra, the mean and most likely maps, and their variances and covariances, can be computed if desired.

Bayesian estimation of the low-energy constants up to fourth order in the nucleon-nucleon sector of chiral effective field theory

We use Bayesian methods and Hamiltonian Monte Carlo (HMC) sampling to infer the posterior probability density function (PDF) for the low-energy constants (LECs) up to next-to-next-to-next-to-leading order (N3LO) in a chiral effective field theory ($\ensuremath{\chi}\mathrm{EFT}$) description of the nucleon-nucleon interaction. In a first step, we condition the inference on neutron-proton and proton-proton scattering data and account for uncorrelated $\ensuremath{\chi}\mathrm{EFT}$ truncation errors. We demonstrate how to successfully sample the 31-dimensional space of LECs at N3LO using a revised HMC inference protocol. In a second step we extend the analysis by means of importance sampling and an empirical determination of the neutron-neutron scattering length to infer the posterior PDF for the leading charge-dependent contact LEC in the $^{1}S_{0}$ neutron-neutron interaction channel. While doing so we account for the $\ensuremath{\chi}\mathrm{EFT}$ truncation error via a conjugate prior. We use the resulting posterior PDF to sample the posterior predictive distributions for the effective range parameters in the $^{1}S_{0}$ wave as well as the strengths of charge-symmetry breaking and charge-independence breaking. We conclude that empirical point-estimate results of isospin breaking in the $^{1}S_{0}$ channel are consistent with the PDFs obtained in our Bayesian analysis and that, when accounting for $\ensuremath{\chi}\mathrm{EFT}$ truncation errors, one must go to next-to-next-to-leading order to confidently detect isospin breaking effects.

Bayesian Seismic Source Inversion With a 3‐D Earth Model of the Japanese Islands

AbstractWe present probabilistic centroid‐moment tensor solutions inferred from the combination of Hamiltonian Monte Carlo sampling and a 3‐D full‐waveform inversion Earth model of the Japanese islands. While the former provides complete posterior probability densities, the latter allows us to exploit waveform data with periods as low as 15 s. For the computation of Green's functions, we employ spectral‐element simulations through the radially anisotropic and visco‐elastic model, leading to substantial improvements of data fit compared to layered models. Focusing on 13 Mw 4.8–Mw 5.3 offshore earthquakes with a significant non‐double‐couple (non‐DC) component, we simultaneously infer the centroid location, time and moment tensor without any a priori constraints on the faulting mechanism. Furthermore, we perform the inversions across several period bands, varying the minimum period between 15 and 50 s. Accounting for 3‐D Earth structure at shorter periods can increase the double‐couple (DC) component of an event, compared to the GCMT solution, by tens of percent. This suggests that non‐DC events in the GCMT catalog may result from unmodeled Earth structure and the related limitation to longer‐period data. We also observe that significant changes in source parameters, and the DC component in particular, may be related to only small waveform changes, thereby accentuating the importance of a reliable Earth model. Posterior probability density distributions become increasingly multimodal for shorter‐period data that provide tighter constraints on source parameters. This implies, in our specific case, that stochastic approaches to the source inversion problem are required for periods below ∼20 s to avoid trapping in local minima.

Field-based physical inference from peculiar velocity tracers

ABSTRACT We present a proof-of-concept Bayesian hierarchical modelling approach to reconstruct the initial cosmic matter density field constrained by peculiar velocity observations. Using a model for the gravitational evolution of dark matter to connect the initial conditions to late-time observations, it reconstructs the late-time density and velocity fields as natural byproducts. We implement this field-based physical inference approach by adapting the Bayesian Origin Reconstruction from Galaxies ($\small {\rm BORG}$) algorithm, which explores the high-dimensional posterior through the use of Hamiltonian Monte Carlo sampling. We test the self-consistency of the method using random sets of tracers, and assess its accuracy in a more complex scenario where peculiar velocity tracers are mock haloes drawn from $\small {\rm GADGET2}$ N-body simulations. We find that our framework self-consistently infers the initial conditions, density and velocity fields, and shows some robustness to model mis-specification. Compared with the approach of constrained Gaussian random fields/Wiener filtering, the hierarchical model produces more accurate final density and velocity field reconstructions. It also allows us to constrain the initial conditions by peculiar velocity observations, complementing in this aspect other field-based approaches based on alternative cosmological observables such as galaxy clustering or weak lensing.