Search code examples
numpyscipydistributionscipy.stats

Fitting custom distribution scipy.stats gives overflow


I am trying to fit a generalised error distribution to some data that I have. The form of the distribution is given as enter image description here

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)


Solution

  • 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