4
$\begingroup$

I would like to compute the expectation value $\langle O \rangle = \sum_x P(x) O(x)$ of some random variable over an extremely large sample space that I cannot simply exhaustively go through. Usually I would use Metropolis Monte Carlo for this, but for this particular example I already know the exact normalized probability $P(x)$ with which each sample $x$ occurs.

How can I leverage this? Does this simply mean I do not need to wait for equilibration of my Markov Chain? I feel like knowing the normalization should allow me to use a much more efficient algorithm/method here. Ideally I would like to directly generate samples $x_i$ according to the distribution of $P(x)$ and then compute $1/N \sum_{i=1}^N O(x_i)$ but I do not see an easy way of doing this, apart from MCMC which does not utilise the full information I have available.

Thank you!

$\endgroup$
1
  • $\begingroup$ The advantage of knowing the normalised $P(\cdot)$ is that you can compute how much of the mass has been explored by the Markov chain so far and hence decide on a stopping rule. $\endgroup$
    – Xi'an
    May 4 at 4:33

1 Answer 1

3
$\begingroup$

Another way of looking at the issue of approximating$$\mathfrak I = \sum_{x\in\mathfrak X} p(x)O(x)$$by stochastic techniques is to aim at adding primarily large values of $p(x)O(x)$. Assuming no information is available about the [location or values of the] largest probabilities over $\mathfrak X$, one could consider a self-avoiding Markov chain by never returning to entries in $\mathfrak X$ already visited and move from $X_t$ to $X_{t+1}$ by choosing among neighbouring entries with probabilities proportional to $p(x)O(x)$ or $\exp\{\alpha p(x)O(x)\}$. With the added perk that since $p(x)O(x)$ is computed for these neighbours, they can all be added to the approximation of $\mathfrak I$ and excluded from future steps. The algorithm could stop when the probability of the visited values is close enough to one.

Here is a toy illustration

![enter image description here

where $\mathfrak X$ is an $N\times N$ grid, $O(\cdot)$ is a discretised Normal density (represented by the level set on the above picture) and $\mathfrak X$ is explored by a random exploration until the accumulated probability mass is $0.999$. The white dots in the above picture correspond to the points $x$ that have not been explored.

Here is the attached R code:

#preliminaries
N=5e2;N2=N*N
Y=matrix(0,N,N)#visited points
P=exp(matrix(rnorm(N2),N))
P=P/sum(P)#probability mass function
O=matrix(dnorm((1:N2)%%N,mean=N/2,sd=N/3)+
         dnorm((1:N2)%/%N,mean=N/2,sd=N/3),N)#function O(x)
#stochastic exploration
m=sample(1:N2,1)#current point
t=P[m]*O[m]#targeted expectation
p=P[m]#visited mass
Y[m]=1
while(p<.999){
  b=sample(which(!Y),4)#neighbours
  Y[b]=1
  p=p+sum(P[b])
  t=t+sum(P[b]*O[b])
  m=sample(b,1,prob=P[b]*O[b])}#next value
$\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.