Whenever the invariant stationary density of metastable dynamical systems decomposes into almost invariant partial densities, its computation as eigenvector of some transition probability matrix is an ill-conditioned problem. In order to avoid this computational difficulty, we suggest applying an aggregation/disaggregation method which addresses only well-conditioned subproblems aud thus results in a stable algorithm. In contrast to existing methods, the aggregation step is done via a sampling algorithm which covers only small patches of the sampling space. Finally, the theoretical analysis is illustrated by two biomolecular examples.