Hi everyone, I have a program in Fortran and it relies on the use of Gaussian numbers, which are generated within the program. Suppose I now wish to sample numbers from a normal distribution with skewness and kurtosis, is there a way I can manipulate the random numbers generated to do so?

I have been pointed to the Cornish Fisher type expansion:
Z'=Z+(g1/6)(Z^2-1)+(g2/24)(Z^3-3Z)-(g1^2/36)(2Z^3-5Z)

where Z ~ N(0,1). I also know the range of valid values for g1 and g2 which make Z' a valid distribution. However when I implement this method the results are very unstable even with 1000000 simulated numbers. I was wondering if there was any other ways to accomplish what i need?

I also have been pointed towards using an expansion of the form:

Z'=a+bZ+cZ^2+dZ^3

I then calculated up to the 4th moment of this and solved for a,b,c,d such that the mean was 0 and variance was 1 as required and then various amounts of skewness(-2..2) and kurtosis(0..11). I then used these values in the Fortran program but again results seemed unstable, although they did seem relatively stable for the negative skewness values.

Any points, further reading, or alternative methods would be much much appreciated.

Jamie.

Decide what distribution you actually want and generate the numbers directly.