4
$\begingroup$

I'm a phd student in Geography, i need some help (or good ressources) to understand why and when i need to use PIT (Probability integral transform) in my validation program for simulation.

I explain the context : Actually, i'm running a lot of stochastic simulations (more than 100000) because i use genetic algorithm to optimize parameters values in my model.

In this model, we work on system of settlements in old times, and we try to obtain the best hierarchy of settlements. To obtain this result we construct a specific function fitness with multi goal objective.

One of this objective is really strict, i take all of my population at the end of the simulation and i try to test the possible lognormality of this rank size distribution => result is ok, or not ok...

Rank size distribution

Actually, i use a classic kolmogorov test to validate/invalidate the lognormality of my city population data, but ouch, this distribution is "unknown" at the end of simulation, so i use this step to resolve this problem :

1 - I try to estimate the standard deviation and the mean of my "unknow distribution" with MLE method (http://en.wikipedia.org/wiki/Maximum_likelihood)

2 - With std and mean for paramter, i generate a new lognormal distribution, named "generated theorical distribution"

3 - I make a uniform transformation on my data, using the "Probability integral transform" or PIT method. But here i don't understand this step, why i cannot use directly my data, perhaps you can enligth me on this point ? ...

4 - I compare the "generated theorical distribution" with the sorted simulation data after the uniform transformation (step3)

5 - After that i extract p-value, and D value to validate / invalidate the lognormality of my simulation data.

Now, my question on this method :

  • I'm not statistician,so perhaphs these step are false ... in this case, can you explain me the step with the problem ?
  • Also, i don't understand the step 3, and the necessity of the "Probability integral transform" of my simulation data here, can you enlight me on this point ?
  • Do you think i can use this methodology, and ks.test with only 100 city in my model, so only 100 values in the sample used in lognormality test ?

Update 1

The code which run the PITS :

  def evalKs(sample: Array[Double],dist:LognormalDist)=
  {
    data = new DoubleArrayList(sample)
    v = data.elements
    n = data.size
    val dataUnif = unifTransform (data, dist)
    dataUnif.quickSortFromTo (0, dataUnif.size - 1)
    ret = kolmogorovSmirnov (dataUnif)  
  }

The code for unifTransform in SSJ library GOFStat ( http://www.iro.umontreal.ca/~simardr/ssj/doc/html/umontreal/iro/lecuyer/gof/GofStat.html ) :

   public static DoubleArrayList unifTransform (DoubleArrayList data,
                                        DiscreteDistribution dist) {
       double[] v = data.elements();
       int n = data.size();

       double[] u = new double[n];
       for (int i = 0; i < n; i++)
          u[i] = dist.cdf ((int)v[i]);
       return new DoubleArrayList (u);
   }

In the class LognormalDist (by default mu = 0, sigma = 1) :

   public double cdf (double x) {
      return cdf (mu, sigma, x);
   }

   public static double cdf (double mu, double sigma, double x) {
      if (sigma <= 0.0)
        throw new IllegalArgumentException ("sigma  <= 0");
      if (x <= 0.0)
         return 0.0;
      return NormalDist.cdf01 ((Math.log (x) - mu)/sigma);
   }

In the class NormalDist

   public static double cdf01 (double x) {
   /*
    * Returns P[X < x] for the normal distribution.
    * As in J. L. Schonfelder, Math. of Computation, Vol. 32,
    * pp 1232--1240, (1978).
    */

      double t, r;

      if (x <= -XBIG)
         return 0.0;
      if (x >= XBIG)
         return 1.0;

      x = -x/Num.RAC2;
      if (x < 0) {
         x = -x;
         t = (x - 3.75) / (x + 3.75);
         r = 1.0 - 0.5 * Math.exp ( -x * x) * Num.evalCheby (NORMAL2_A, COEFFMAX, t);
      } else {
         t = (x - 3.75) / (x + 3.75);
         r = 0.5 * Math.exp ( -x * x) * Num.evalCheby (NORMAL2_A, COEFFMAX, t);
      }
      return r;
   }

Thanks a lot for your advice, critics on this students work in progress :) SR

$\endgroup$
  • $\begingroup$ Reyman, could you post your code and a sample of the way you calculate the PITS? $\endgroup$ – AndresT Feb 13 '12 at 22:40
  • $\begingroup$ I update my post with code :) $\endgroup$ – reyman64 Feb 14 '12 at 16:03

Your Answer

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

Browse other questions tagged or ask your own question.