CODEWITHSUNDEEP
pythonnumpyscipyprobability-distribution
baccandr
baccandr

Reputation: 1130

How to fit a distribution defined with scipy.stats.rv_continuous?

I would like to fit data with a combination of distributions in python and the most logical way it seems to be via scipy.stats.rv_continuous. I was able to define a new distribution using this class and to fit some artificial data, however the fit produces 2 variables more than the free parameters of the distribution and I don't understand how to interpret these. In addition, the fit is very slow so any suggestion on how to speed it up would be highly appreciated.

Here a minimum reproducible example (for the sake of this question I will be using the combination of a normal and a lognormal distributions):

import numpy as np
import scipy.stats as stats

# Create the new distribution combining a normal and lognormal distr
def lognorm(x,s,loc,scale):
    return(stats.lognorm.pdf(x, s = s, loc = loc, scale = scale))
def norm(x,loc,scale):
    return(stats.norm.pdf(x, loc = loc, scale = scale))

class combo_dist_gen(stats.rv_continuous):
    "Gaussian and lognormal combination"
    def _pdf(self, x,  s1, loc1, scale1, loc2, scale2):
        return (lognorm(x, s1, loc1, scale1) + norm(x, loc2, scale2))

combo_dist = combo_dist_gen(name='combo_dist')

# Generate some artificial data
gen_data = np.append(stats.norm.rvs(loc=0.2, scale=0.1, size=5000),\
    stats.lognorm.rvs(size=5000, s=0.1, loc=0.2, scale=0.5))

# Fit the data with the new distribution
# I provide initial values not too far from the original distribution
Fit_results = combo_dist.fit(gen_data, 0.15, 0.15, 0.6, 0.25, 0.05)

A part from being very slow the fit seems to work, however it returns 7 variable while the original distribution only has 5 free parameters:

print(Fit_results)
(0.0608036989522803, 0.030858042734341062, 0.9475658421131599, 0.4083398045761335, 0.11227588564167855, -0.15941656336149485, 0.8806248445561231)

I don't understand what these 2 additional variables are and how they enter into the definition of the distribution.

If I generate a new pdf using the fit results I can reproduce well the original distribution but only using all the 7 variables:

xvals = np.linspace(-1,3, 1000)
gen_data_pdf = (lognorm(xvals,0.1, 0.2, 0.5)+norm(x, 0.2,0.1))
ydata1 = combo_dist.pdf(xvals,*Fit_results)
ydata2 = combo_dist.pdf(xvals,*Fit_results[:5])

plt.figure()
plt.plot(xvals, gen_data_pdf, label = 'Original distribution')
plt.plot(xvals, ydata1, label = 'Fitted distribution, all parameters')
plt.plot(xvals, ydata2, label = 'Fitted distribution, only first 5 parameters')

plt.legend()

Here the result

p.s.1 The official documentation is a bit obscure to me and doesn't seem to provide any useful example. Here on SO there are a few answers providing some explanations (like here and here) but none of them seem to address my issue.

p.s.2 I am aware that the pdf of the combined distribution is not normalized to 1. In my original implementation I was dividing the pdf by 2 but for some reason with the additional division the fit didn't work (RuntimeError, no convergence)

Upvotes: 2

Views: 1362

Answers (1)

danifos
danifos

Reputation: 61

The 2 variables are the loc and scale parameters to shift and scale the distribution according to the documentation. Just fix the values by:

Fit_results = combo_dist.fit(gen_data, 0.15, 0.15, 0.6, 0.25, 0.05,
                             floc=0, fscale=1)

Upvotes: 2

Related Questions