[Bug-gsl] levy skew alpha-stable distribution documentation mis-plot?

[dorkusmonkey]

  On p. 218 for my copy of the gsl reference, chapter 19.13 "The Levy
skew alpha-Stable Distribution" I noticed that the plot of the
probability distribution function, p(x), (for alpha=1.0, beta=1.0,
scale=1.0) has a maximum at around y=0.06 or so.  I have estimated the
pdf myself and have found the peak for the same parameters to be
nearer to 0.3.  Are you sure that the plot of the pdf (p(x)) is
correct?

   below is a small program that estimates the integral and spits out
the pdf in the range of [-10,10] in gnuplot friendly output, along
with the cdf in comments.  I hope that I have interpreted everything
correctly and have done the integration and plot correctly, but I
would not put it past me to mess it up, so I would greatly appreciate
any clarifications of misconceptions/fallacies I might have.

Thank you,
dorkusmonkey


----------------------------------------
/* gcc -lm levy_pdf.c -o levy_pdf */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define sign(x) ( ((x)<0.0) ? -1.0 : 1.0 )

/* envelope (kernel?) */
double env(double b, double x) {
  if ((x<-b) || (x>b)) return(0.0);
  x = sign(x)*x;
  return( 1.0 - (x/b) );
}

/* special case when alpha=1 */
double levy_pdf_a1(double x, double gamma, double beta,
    double b, double dx) {
  int i, j, k, n;
  double t, y, s=0.0, ct_a;
  for (n=0, t=-b; t<=b; t=-b + ((double)(++n)*dx)) {
    ct_a = fabs(gamma*t);
    s += env(b, t) * (exp(-ct_a)
      *  cos((t*x) - (beta* (-2.0/M_PI)*log(fabs(t))*ct_a*sign(t))) )*dx;
  }
  return(s / (2.0*M_PI));
}


double levy_pdf(double x, double alpha, double gamma, double beta,
    double b, double dx) {
  int i, j, k, n;
  double t, y, s=0.0, ct_a;
  if (fabs(1.0-alpha)<0.001) return(levy_pdf_a1(x, gamma, beta, b, dx));
  for (n=0, t=-b; t<=b; t=-b + ((double)(++n)*dx)) {
    ct_a = pow(fabs(gamma*t), alpha);
    s += (exp(-ct_a)
      *  cos((t*x) - (beta*tan(M_PI*alpha/2.0)*ct_a*sign(t))) )*dx;
  }
  return(s / (2.0*M_PI));
}

int main(int argc, char **argv) {
  double x, y, dx, lx, rx;
  int i, j, k, n;

  double alpha, gamma, beta;
  double s=0.0, t;


  alpha = 1.0;
  gamma = 1.0;
  beta = 1.0;

  dx=0.01;

  lx = -10.0;  rx = 10.0;

  for (n=0, x=lx; x<=rx; x = lx + ((double)(++n)*dx)) {

    t = levy_pdf(x, alpha, gamma, beta, 100.0, 0.001);
    printf("%g %g\n", x, t);

    s += t*dx;

    printf("# %g %g\n", x, s);

  }

}
-----------------------------------------