Search code examples
mathstatisticsgeometryfortrancomputational-geometry

Drawing random values from a Fisher Distribution


In my research, I am generating discrete planes that are intended to represent fractures in rock. The orientation of a fracture plane is specified by its dip and dip direction. Knowing this, I also know the components of the normal vector for each plane.

So far, I have been drawing dip and dip direction independently from normal distributions. This is OK, but I would like to add the ability to draw from the Fisher distribution.

The fisher distribution is described HERE

Basically, I want to be able to specify an average dip and dip direction (or a mean vector) and a "fisher constant" or dispersion factor, k, and draw values randomly from that orientation distribution.

Additional Info: It seems like the "Von Mises-Fisher distribution" is either the same as what I've been calling the "Fisher distribution" or is somehow related. Some info on the Von Mises-Fisher distribution:

As you can see, I've done some looking into this, but I admit that I don't fully understand the mathematics. I feel like I'm close, but am not quite getting it... Any help is much appreciated!

If it helps, my programming is in FORTRAN.


Solution

  • The algorithm is on page 59 of "Statistical analysis of spherical data" by N. I. Fisher, T. Lewis and B. J. J. Embleton. I highly recommend that book -- it will help you understand the mathematics.

    The following will produce random Fisher distribution locations centered on the North pole. If you want them randomly centered, then you produce additional uniform random locations on the sphere and rotate these locations to be centered on those locations. If you are not sure of those steps, consult the aforementioned book. This Fortran code fragment uses a random number generator that produces uniform deviates from 0 to 1.

      lambda = exp (-2.0 * kappa)
      term1 = get_uniform_random () * (1.0 - lambda) + lambda
      CoLat = 2.0 * asin ( sqrt ( -log (term1) / (2.0 * kappa) ) )
      Long = 2.0 * PI * get_uniform_random ()