1

I'm trying to replicate the answer from R's gamlss.dist exGAUS in Python's scipy.stats.exponnorm.

The following in R returns -0.4003516:

library(gamlss)
qexGAUS(0.5, mu = -1., sigma = .6, nu = .7, lower.tail = TRUE, log.p = FALSE)

However, the following in scipy returns -0.6153701019503552:

from scipy.stats import exponnorm
exponnorm.ppf(q=0.5, K=.7, loc=-1, scale=.6)

How do I configure the scipy inputs to get the same answer as I get in R? Based on Why does gamlss give incorrect estimates of ex-gaussian distribution parameters?, I think it is related to log transforms of the inputs, but I can't seem to get a match with various combinations of log/antilog inputs.

a11
  • 3,122
  • 4
  • 27
  • 66

1 Answers1

2

As you surmised, the parameterizations for the two functions are different.

Both functions have three parameters. In R the parameters are mu (mean of normal), sigma (standard deviation of normal), and nu (mean of the exponential). In Python the parameters are mu (mean of the normal), sigma (standard deviation of the normal), and K (shape of the exponential).

The Python help documentation tells us that in order to translate from one to the other we need to use K = 1/(sigma*lambda).

In this formulation lambda is the rate or inverse scale of the exponential meaning lambda = 1/nu. Therefore, we need to use K = 1/(sigma*(1/nu)) to translate from R to Python:

from scipy.stats import exponnorm
print(exponnorm.ppf(q=0.5, K=1/(.6*(1/.7)), loc=-1, scale=.6))

The above produces -0.4003398910762068, where the difference at the 5th decimal place is likely due to differing estimation routines.

the-mad-statter
  • 5,650
  • 1
  • 10
  • 20