2
$\begingroup$

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.

enter image description here

$\endgroup$
  • $\begingroup$ 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 integrate $p/q$. $\endgroup$ – Kirill May 29 '15 at 1:09
  • $\begingroup$ thanks. To the second part - do you think that this approach is somehow usefull for solving problems like the two above? I mean in this case whole problem is actually to compute the total integral, which is basically the normalization constant itself. So If I'm correct If I would know some estimate of the total integral from other calculation Metropolis would not improve it any further. ( ? ) $\endgroup$ – Prokop Hapala May 29 '15 at 8:05
  • $\begingroup$ I expanded my comment into an answer. $\endgroup$ – Kirill May 29 '15 at 18:03
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I like the idea about sampling more densly in the centre of plates and than renormalizing by that proposal distribution. This is probably what I was searching for. I was obsessed with miss-guiding idea that preferential sampling = Metropolis algorithm. $\endgroup$ – Prokop Hapala May 29 '15 at 21:08
  • $\begingroup$ ad solid angle - In principle I know how to compute the solid angle from point either by Oosterom and Strackee from wiki or by an other Monte Carlo. But I find it to laborious ( both for me as programer, as well as for number CPU operations required - considering that is should be computed for every $i$ ). Also I wanted to learn some general tool for solution of similar problems, which would not require problem-specific optimization. Metropolis-Hasting did sound like such general tool. But perhaps for other class of problems. $\endgroup$ – Prokop Hapala May 29 '15 at 21:11
  • $\begingroup$ @ProkopHapala You can precompute it (it only depends on the first point's position relative to the second plate), using interpolation techniques, making it cheaper. $\endgroup$ – Kirill May 29 '15 at 23:26

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.