Estimate Ratio of Normalizing Constants from two datasets

Question

Suppose I have a non-negative function $f:\mathbb{R}^N \to [0, +\infty)$ that defines two different (unnormalized) probability densities on two separate subsets $A, B \subset \mathbb{R}^N$ with $A \cap B = \emptyset$.

$$ \widetilde{p}(x) = \begin{cases} f(x) & \text{if } x\in A \\ 0 & \text{if } x \notin A\end{cases} \\ \widetilde{q}(x) = \begin{cases} f(x) & \text{if } x\in B \\ 0 & \text{if } x\notin B\end{cases} $$ Now suppose a set of samples from each $$ x_1^p,\ldots, x_N^p \sim \widetilde{p}(x) \\ x_1^q, \ldots, x_N^q \sim \widetilde{q}(x) $$ How can I use these samples to approximate the ratio of normalizing constants $$ \frac{Z_q}{Z_p} = \frac{\displaystyle \int_B \widetilde{q}(x)dx}{\displaystyle \int_A \widetilde{p}(x)dx} $$

Answer 1

[Warning: the following proposals, except one, assume $f$ can be evaluated at arbitrary values. If the only inputs are the sets $A$ and $B$, and the samples $\mathbf x$ and $\mathbf y$, the problem has no solution since any combination $\alpha f_A(x)+(1-\alpha) f_B(x)$ could produce exactly the same samples $\mathbf x$ and $\mathbf y$.]

If an expectation under $f$ [assumed to be a density, i.e. normalised] is known, as in control variate settings, $$\int h(x) f(x)\,\text dx=\mathfrak h>0$$ solving$$\alpha\int_A h(x) \frac{\tilde p(x)}{Z_p} \,\text dx+(1-\alpha) \int_B h(x) \frac{\tilde q(x)}{Z_q} \,\text dx=\mathfrak h\tag{1}$$returns $$\alpha=Z_p$$unless both integrals are identical. Replacing (1) with its Monte Carlo version $$\sum_{i=1}^N\{\hat\alpha h(x_i)+(1-\hat\alpha)h(y_i)\}=n\mathfrak h$$ returns an estimator of $Z_p$. In the case $f$ is not normalised, two such functions $h$ would be required.

An illustration in the Normal $\mathcal N(0,1)$ case when $A=(-\infty,1)$, $p_A=0.841$, $h(x)=x$, and $\mathfrak h=0$:

for(t in al<-1:1e2){
  x=(x<-rnorm(4e5/p))[x<1][1:1e5]
  y=(y<-rnorm(4e5/(1-p)))[y>1][1:1e5]
  al[t]=1/(1-mean(x)/mean(y))}

resulting in a concentrated approximation of $p_A$:

Another (crude) approach would be to estimate (by a kernel density estimator) the densities on both samples and, assuming continuity of the density $f$ at the boundary, make them meet at this boundary (between $A$ and $B$). As in this example

p=pnorm(1) #A is (-oo,1) and B (1,+oo) and f is dnorm
wh=N<-1:1e5
x=(x<-rnorm(4*N/p))[x<1][1:N] #sample from p over A
y=(y<-rnorm(4*N/(1-p)))[y>1][1:N] #sample from q over B
a=density(x,ker="gaussian",to=1) #KDE of p
hatfa=approxfun(a$x,a$y)
b=density(y,ker="gaussian",from=1) #KDE of q
hatfa=approxfun(a$x,a$y)

leading to

> hatfb(1)/hatfa(1)
[1] 4.917869

as the estimate of $p_A/p_B$. Repeated simulations show however that the estimator is biased downwards.

As an alternative, once one is using a kernel approximation, one can implement a full (reversible jump) MCMC version for simulating moves between $A$ and $B$, based on the non-parametric density estimator as a proposal:

mm=wh=1:N;wh[1]=x[1] #mm model indicator, 1 stands for A
for (t in 2:N){#Metropolis-Hastings steps
  mm[t]=mm[t-1];wh[t]=wh[t-1]
  if(mm[t]){#propose to move to B
    prop=rnorm(1,sample(y,1),b$bw)
    while(prop<1)prop=rnorm(1,sample(y,1),b$bw)#constrained to B
    if(dnorm(wh[t])*hatfb(prop)*runif(1)<dnorm(prop)*hatfa(wh[t])){
      wh[t]=prop;mm[t]=0}
  }else{#propose to move to A
    prop=rnorm(1,sample(x,1),a$bw)
    while(prop>1)prop=rnorm(1,sample(x,1),a$bw)#constrained to A
    if(dnorm(wh[t])*hatfa(prop)*runif(1)<dnorm(prop)*hatfb(wh[t])){
      wh[t]=prop;mm[t]=1}}
}

resulting in an estimate of the ratio

> mean(mm)/mean(!mm)
[1] 5.127451