Python: comment définir des distributions personnalisées?

Question

Je veux créer une distribution personnalisée basée sur une loi tronquée de Levy, qui lit

p(r) = (r + r0)**(-beta)*exp(-r/k).

I Il définit de la manière suivante:

import numpy as np 
import scipy.stats as st 

class LevyPDF(st.rv_continuous): 
    def _pdf(self,r): 
     r0 = 100 
     k = 1500 
     beta = 1.6 
     return (r + r0)**(-beta)*np.exp(-r/k) 

Supposons que je veux trouver la répartition des distances entre r = 0 et r = 50km. Puis:

nmin = 0 
nmax = 50 
my_cv = LevyPDF(a=nmin, b=nmax, name='LevyPDF') 
x = np.linspace(nmin, nmax, (nmax-nmin)*2) 

Je ne comprends pas pourquoi:

sum(my_cv.cdf(x)) = 2.22 

au lieu de 1.

Alors, comment puis-je définir un histogramme de N = 2000000 distances aléatoires sur la base de la distribution que je définissais?

Answer 1

En utilisant votre exemple minimal (légèrement adapté):

import scipy.stats as st 
import numpy as np 
import matplotlib.pyplot as plt 

class LevyPDF(st.rv_continuous): 
    def _pdf(self,r): 
     r0 = 100 
     k = 1500 
     beta = 1.6 
     return (r + r0)**(-beta)*np.exp(-r/k) 

nmin = 0 
nmax = 50 
my_cv = LevyPDF(a=nmin, b=nmax, name='LevyPDF') 

Pour échantillon à partir de votre variable aléatoire, utilisez la méthode rvs() de rv_continuous classe:

N = 50000 
X = my_cv.rvs(size=N, random_state=1) 

retourne un tableau de taille (N,) avec aléatoire variates échantillonnées à partir de votre distribution. Utilisez l'option random_state pour geler votre exemple et rendre votre script reproductible (il définit une graine aléatoire pour votre échantillonnage).

note que N augmente doucement, le temps de calcul augmente de façon drastique.

Pour tracer l'histogramme, utilisez la bibliothèque matplotlib, voir hist:

fig, axe = plt.subplots() 
n, bins, patches = axe.hist(X, 50, normed=1, facecolor='green', alpha=0.75) 
plt.show(axe) 

Bellow un exemple d'échantillonnage de Chi carré avec 40 degrés de liberté:

from scipy import stats 
import numpy as np 
import matplotlib.pyplot as plt 
rv = stats.chi2(40) 
N = 200000 
X = rv.rvs(size=N, random_state=1) 
fig, axe = plt.subplots() 
n, bins, patches = axe.hist(X, 50, normed=1, facecolor='green', alpha=0.75) 
plt.show(axe) 

Il conduit: