# Estimate Ratio of Normalizing Constants from two datasets

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\int_B \widetilde{q}(x)dx}\int_A \widetilde{p}(x)dx}$$

• @Xi'an what an honor! Thank you! I had a look at some of those papers like the exchange algorithm etc but they all use some other proposal density defined on the same space. In my case, however, I have two disjoint supports. Or rather, the supports are not disjoint but one of the densities is always $0$ when the other is positive Apr 27 at 15:48
• You are correct, this is not a straightforward application of the above. Since the simulations are constrained to $A$ and $B$ respectively, they bring no information about the relative weights of $A$ and $B$ under $f$. This reminds me of the difficulty of computing the Bayes factor when given only samples from each posterior. Apr 28 at 5:11
• Do you have an expression for $f$, or just a black-box sampler for $\tilde p$ and $\tilde q$? Would you be able to sample from the density that is proportional to $f$ over another subset of $\mathbb{R}^N$? Apr 28 at 17:25
• @RobinRyder unfortunately, I don't have an expression for $f$. In fact, I would like to learn about $Z_p$ and $Z_q$ so that I can then learn about $f$! Apr 29 at 10:48
• So evaluating or estimating $f(a)/f(b)$ with $a\in A, b\in B$ would not be possible? Apr 29 at 10:52

[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)


> hatfb(1)/hatfa(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=x #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)
 5.127451