I am trying to fit a generalised error distribution to some data that I have. The form of the distribution is given as
I have tried the following implementation
import numpy as np
import scipy.stats as st
from scipy.special import gamma
class ged(st.rv_continuous):
def _pdf(self, x, mu, sigma, kappa):
term1 = gamma(3*kappa)/gamma(kappa)
exponent = (term1 * ((x - mu)/sigma)**2)**(1/(2*kappa))
term2 = np.exp(-exponent)
term3 = 2*sigma*gamma(kappa+1)
fx = term1**0.5 * term2/term3
return fx
ged_inst = ged(name='ged')
data = np.random.normal(size=1000)
ged_inst.fit(data, 0, 0.01, 1)
However this gives
OverflowError: (34, 'Numerical result out of range')
How do I correctly implement this distribution? I am trying to fit to real data (not the toy normal data generated in the question)
As posted in the comments, to get this working I needed to override the default _argcheck
function. The following works:
class ged(st.rv_continuous):
def _pdf(self, x, mu, sigma, kappa):
term1 = gamma(3*kappa)/gamma(kappa)
exponent = (term1 * ((x - mu)/sigma)**2)**(1/(2*kappa))
term2 = np.exp(-exponent)
term3 = 2*sigma*gamma(kappa+1)
fx = term1**0.5 * term2/term3
return fx
def _argcheck(self, mu, sigma, kappa):
s = sigma > 0
k = kappa < 1
return s and k