7
$\begingroup$

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} $$

$\endgroup$
5
  • 1
    $\begingroup$ @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 $\endgroup$ Apr 27 at 15:48
  • 1
    $\begingroup$ 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. $\endgroup$
    – Xi'an
    Apr 28 at 5:11
  • $\begingroup$ 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$? $\endgroup$ Apr 28 at 17:25
  • $\begingroup$ @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$! $\endgroup$ Apr 29 at 10:48
  • 1
    $\begingroup$ So evaluating or estimating $f(a)/f(b)$ with $a\in A, b\in B$ would not be possible? $\endgroup$ Apr 29 at 10:52
6
$\begingroup$

[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$:

enter image description here

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.