Metropolis Monte Carlo integration of Area with unknown normalization

Question

I probably miss something very basic. I don't see how to use Metropolis–Hastings algorithm for computation of integrals, if I don't know the volume of accessible phase space (i.e. proper normalization by partition function).

Consider this two problems ( which are very much connected):

1) Number of possible states under some constrains

I have some phase-space with some areas possible ( probability $p(X) = 1$ ) and some impossible ( probability $p(X) = 0$ ). I don't know the volume of neither areas, but I can easily check if a point is in possible or inpossible area.

Because I know that area of impossible part of phase-space is much bigger than area of possible part, I wan't to use Metropolis–Hastings algorithm to preferentially explore the possible configurations. I can easily implement the updates according to Metropolis–Hastings to produce Markov Chain which does explore the possible part of phase-space, however, I don't see how to use it to actually compute the possible area by summing over such Markov chain. I still don't know proper normalization.

There is a javascript visualization of the problem in 2D http://www.openprocessing.org/sketch/200710 the most relevant part of code is there:

// in this example probability of configuration is p=0 
// if the x,y point overlap with some sphere, otherwise it is p=1
float getUnnormalizedProbability( float x, float y ){
  float prob = 1.0;
  for(int i=0; i<nsph; i++){  
    float dx = x - sphx[i];
    float dy = y - sphy[i];
    float r2 = dx*dx + dy*dy;
    if( r2 < sphR2[i] ){
      prob = 0.0;
      break;
    }
  }
  return prob;
}

void MMCstep(){
  float newx = x + random(-stepsz,stepsz);
  float newy = y + random(-stepsz,stepsz);
  // Periodic boundary condition
  if ( newx<0      ) newx = boxsz + newx;
  if ( newx>boxsz  ) newx = newx  - boxsz;
  if ( newy<0      ) newy = boxsz + newy;
  if ( newy>boxsz  ) newy = newy  - boxsz;
  float newp = getUnnormalizedProbability( newx, newy );
  // Metropolis-Hastings condition  alfa =  min( 1, pnew/p );
  if( newp >= p ){
    move( newx, newy, newp );
  }else{
    if( (newp/p) > random(1.0) ){
      move( newx, newy, newp );
    }
  }
}

2) Neutron scattering on finite plates

I have source of radiation $S$ and detector $D$ (of some size) and several thin plates of neutron scattering material (of some size and with some normalized angular scattering distribution).

I want to compute percentage of neutrons emitted by $S$ which impact on $D$.

It can be computed by path-integration over paths by randomly choosing points $i$ and $j$ on plates. This however require renormalization of the sum by probability that a neutron originating in point $i$ impact second plate at all, because I don't sample situations where neutron does not hit the plate, I don't know this renormalization.

Explicity analytical computation of this probability renormalization (i.e. solid angle under which plat_2 is visible from point $i$ ) is quite laboriously. I don't want to do that.

I wonder if Metropolis–Hastings algorithm can help with this somehow.

Answer 1

I don't think you're missing anything, MCMC is used to sample points from a given distribution, known up to a normalization constant, and to evaluate expectations w.r.t. that distribution (not integrals over the volume).

It's often used when the normalization constant of the distribution is not known, and doesn't require you to know it. On the other hand, if you had two distributions, you could estimate the ratio of their normalizing constants. Since you can pick any MC, you can pick one with a different stationary distribution $q$ with a known normalization constant, then use that to estimate $\mathbb{E}^Q[p/q]$ which would compute the integral $\int p(x)\,dx$ over the domain of $q$.

As to your comment, I can think of a few things to try, but I find it difficult to be definitive here. I think at some point you might just have to compute the probability that a neutron hits plate 2.

Simulating paths that don't hit plate 1 or hit plate 1 but don't hit plate 2 should work: depending on how big the plates are, this might work just fine. This should let you know analytically the appropriate normalization constant. I don't see too much reason to avoid this at all costs: for one thing, whatever your final method is, you could also compute this probability by letting $D$ have infinite size and letting $S$ and plate 1 have point sizes and maybe close to each other, then the probability is mostly the probability of plate 2 being hit from the point plate 1.

If you do sample points that don't lie on plate 2, you don't have to sample them uniformly from all of the plane in which it lies: e.g., you can pick some Gaussian distribution approximately centered on top of plate 2.

Also, if computing the solid angle is really the main thing stopping you from normalizing your functions, then it should be possible to do it numerically. Looking at Wikipedia, there is a formula for solid angle subtended by a triangle with given vertices, so to calculate the solid angle subtended by plate 2 from point $i$, you can just split plate 2 into a suitable number of triangles and sum up individual solid angles for each triangle. This might be the most straightforward solution here.