$\newcommand\th\theta\newcommand\R{\mathbb R}$Suppose that we have a parametric setting, that is, the unknown distribution $P$ belongs to a known parametric family $(P_\th)$ of distributions parameterized by a sufficiently low-dimensional parameter $\th$; this may be the case if $P$ is Gaussian. 
Then you can get (say) a maximum likelihood estimate (MLE) $\hat\th_n$ of $\th$ and estimate $R$ by 
$$\hat R_n:=\int_H p_{\hat\th_n}\,ds,$$
where $p_{\hat\th_n}$ is the density of $P_{\hat\th_n}$. 

The MLE is usually consistent, so that you will have $\hat\th_n\to\th$ in probability (as $n\to\infty$). If now $p_\th$ is  continuous on $\th$ uniformly on compact subsets of $\R^m$, you will get $p_{\hat\th_n}\to p_\th$ uniformly on compact subsets of $\R^m$.  If, moreover, you can control the tails of the densities $p_\th$ so as to have their local uniform integrability (as you would have in the Gaussian case), then you will end up with the desired conclusion $\hat R_n\to R$ in probability. 

--- 

If the setting is nonparametric but the dimension $m$ is small, then you can use nonparametric (say kernel) estimators of the unknown density $p$, instead of the parametric estimators. However, then, naturally, the convergence $\hat R_n\to R$ will be substantially more problematic. 

If the setting is completely nonparametric and the dimension $m$ is not small (say $m\ge10$), then the situation will of course be even worse. Here you will need the sample size $n$ to be at least as big as something like $100^m\ge10^{20}$.