CODEWITHSUNDEEP
pythonscipystatisticsdistribution
develarist
develarist

Reputation: 1375

Johnson Moments distribution in Python

The Johnson Moments distribution, whose algorithm was published in 1976, is implemented in

In Scipy, there is no implementation of it, only the Johnson-SU and -SB distributions, which are not the same as the Johnson Moments distribution. Is there another python library, or how can it be implemented in Python?

Hill, I. D, R. Hill, and R. L. Holder. 1976. Algorithm AS99: Fitting Johnson curves by moments. Applied Statistics 25 (2): 180--189

Upvotes: 1

Views: 899

Answers (1)

LTZJoe
LTZJoe

Reputation: 86

I faced the same problem in a side project I'm currently working, and I found a solution that works quite well.

First off, I suggest you read the Applied Statistics article where the algorithm was first published. You'll need to sign up for a free account, but then you get 100 free articles per month. The reason I suggest this is, based on your description, I get the feeling you may have misunderstood what the algorithm does. You mention the "Johnson Moments distribution" which isn't a thing. The algorithm described in that article (same as you mention in your post) describes a function named JNSN that takes mean, stdDev, skewness and kurtosis as inputs, and returns an estimated Johnson distribution type (Su, Sb, Normal or Exponential) plus the 4 parameters needed for the distribution (γ, δ, β and λ). These four parameters needed are described in the literature about the Johnson distributions and are named gamma, beta, sigma and lam respectively (though in the code he refers to the last 2 as XLAM and XI, and has them in the opposite order in his function signatures).

Given the output "itype", you can choose to instantiate either a normal curve, exponential curve, Johnson-SU or Johnson-SB via SciPy. When you do this, beta and gamma correspond to SciPy's "a" and "b" parameters, whereas sigma and lam correspond to the "location" and "scale" parameters. You can test this by pulling the Mean, StdDev, Skewness and Kurtosis from the SciPy distribution you instantiated and checking to see that they match the inputs you passed to the JNSN function.

Now for the python implementation... There isn't one. The code there is fortran. I tried translating that to Python but had too many bugs in my translation. Also, the speed of the translated code was terrible. So, i abandoned the idea of translating the code and instead used numpy's excellent F2PY module to compile the fortran source code into a machine language python plugin.

Note: http://lib.stat.cmu.edu/apstat has the source code to a bunch of other fortran modules as well.
Final Note: The original source code used 4 byte floating numbers. I updated that directly in the fortran code, then used f2py -h ... command to generate, then tweak the signatures for python.

Upvotes: 2

Related Questions