Dimensionality reduction is a key ingredient of many machine learning algorithms and is paramount to their success. For manifold-valued data, the nonlinear equivalent of the well-known principal component analysis (PCA), called, principal geodesic analysis (PGA) is used quite often. An alternative to PGA that is more general and flexible, called ”Nested Homogeneous Spaces (NHS)” for dimensionality reduction of manifold-valued data was recently introduced. In this paper, we present a novel probabilistic version of the NHS model (PNHS) for dimensionality reduction of high dimensional manifold-valued data in Rieman- nian homogeneous spaces. The PNHS model has several advantages over its deterministic counterpart namely, the NHS model. In particular, the ability to, quantify uncertainty in parameter estimates and tackle missing data. We demonstrate these advantages via real and synthetic data examples.