CODEWITHSUNDEEP
pythonnumpyscipyprobability
Dalek
Dalek

Reputation: 4318

Fitting a probability distribution to the data and finding cumulative distribution function for it

I want to fit an asymmetric probability distribution to my data and I thought an exponentially modified Gaussian distribution can be a good representative for my data. I 

m=array([ 16.25,  16.75,  17.25,  17.75,  18.25,  18.75,  19.25,  19.75,
          20.25,  20.75,  21.25,  21.75,  22.25,  22.75,  23.25,  23.75,
          24.25,  24.75,  25.25,  25.75,  26.25,  26.75,  27.25,  27.75,
          28.25,  28.75,  29.25,  29.75,  30.25,  30.75])
pdf=array([  0.00000000e+00,   2.40818784e-04,   1.38470801e-03,
             1.62552679e-03,   3.07043949e-03,   3.37146297e-03,
             5.47862733e-03,   8.36845274e-03,   1.61348585e-02,
             1.92052980e-02,   2.79951836e-02,   3.97953040e-02,
             4.95484648e-02,   7.09211318e-02,   9.50030102e-02,
             1.40878989e-01,   1.90186635e-01,   2.42022878e-01,
             2.77302830e-01,   2.69054786e-01,   2.40397351e-01,
             1.74593618e-01,   9.16917520e-02,   2.41420831e-02,
             7.22456352e-03,   3.01023480e-04,   0.00000000e+00,
             0.00000000e+00,   0.00000000e+00,   6.02046960e-05])

I would like to use scipy.optimize library and in meanwhile could control the goodness of fit maybe and see it in order to improve the chi-squared by changing the initial conditions for input parameters. I have written the following piece of code:

import scipy.special as sse
from math import *
import numpy as np
import scipy.optimize
#defines the PDF of an exponentially modified Gaussian distribution
fitfunc =lambda p,x: 0.5*p[2]*np.exp(0.5*p[2]*(2*p[0]+p[2]*p[1]*p[1]-2*x))*sse.erfc((p[0]+p[2]*p[1]*p[1]-x)/(np.sqrt(2)*p[1]))
"""Deviations of data from fitted curve"""
errfunc = lambda p, x, y: fitfunc(p, x) - y
#initial values
p0=[24,1,1]
p1, success = scipy.optimize.leastsq(errfunc, p0, args=(pdf, m), maxfev=10000)

Update: well I just chose numpy.exp and the first problem solved but still leastsq doesn't give me reliable outputs, what should I do? In addition, I would like to obtain the CDF for this distribution too.

Upvotes: 1

Views: 921

Answers (1)

Robert Dodier
Robert Dodier

Reputation: 17585

The method of least squares IS NOT the method to use for fitting data to a given pdf.

What you (probably) want is the method of maximum likelihood -- i.e. to maximize p(x | a) where a are the parameters of the distribution and x are the data. Typically one forms the log-likelihood and assumes independence, so log p(x | a) = sum(log(pdf(x[i], a)), i, 1, n).

You need to use a minimization function, giving it log p(x | a) as the function to be minimized, and a as its free parameters.

Upvotes: 4

Related Questions