Почему случайные выборки из моего нестандартного дистрибутива не следуют за pdf

Question

Я создал собственный дистрибутив, используя метод scipy rv_continuous. Я пытаюсь создать распределение энергии электрона, произведенного бета-распадом. Учитывая его pdf:

, который я взял из: http://hyperphysics.phy-astr.gsu.edu/hbase/Nuclear/beta2.html#c1

Я определяю свой дистрибутив:

import numpy as np
from scipy.stats import rv_continuous 
import matplotlib.pyplot as plt 

class beta_decay(rv_continuous):
    def _pdf(self, x):
        return (22.48949986*np.sqrt(x**2 + 2*x*0.511)*((0.6-x)**2)*(x+0.511))

# create distribution from 0 --> Q value = 0.6 
beta = beta_decay(a=0, b= 0.6)

# plot pdf 
x = np.linspace(0,0.6)
plt.plot(x, beta.pdf(x))
plt.show()

# random sample the distribution and plot histogram 
random = beta.rvs(size =100)
plt.hist(random)
plt.show()

Где x = KE, Q = 0,6, C = 22,48 ... (определяется путем интегрирования вышеуказанного выражения между 0 -> Q и установки равной 1 для нормализации), и I не учитывать функцию Ферми F (Z ', KEe) в приведенном выше уравнении.

Когда я строю график в формате pdf, он выглядит правильно:

Однако, когда я пытаюсь извлечь из него случайные выборки, используя .rvs (), значение они берутся в массовом порядке по отношению к RHS, а не ниже пика PDF, как я ожидал:

В конечном счете, мой код должен сэмплировать дистрибутив, чтобы получить KE электрона, выпущенного бета-распадом. Почему моя гистограмма такая неправильная?

Answer 1

Я думаю, что ваш PDF определен неправильно, он не нормализован. После того, как я нормализовал его и сделал правильную гистограмму, он, кажется, работает нормально

Код (Win10 x64, Anaconda Python 3.7)

#%%
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate

from scipy.stats import rv_continuous

def bd(x):
    return (22.48949986*np.sqrt(x**2 + 2*x*0.511)*((0.6-x)**2)*(x+0.511))

a = 0.0
b = 0.6

norm = integrate.quad(bd, a, b) # normalization integral
print(norm)

class beta_decay(rv_continuous):
    def _pdf(self, x):
        return bd(x)/norm[0]

# create Q distribution in the [0...0.6] interval
beta = beta_decay(a = a, b = b)

# plot pdf
x = np.linspace(a, b)
plt.plot(x, beta.pdf(x))
plt.show()

# sample from pdf
r = beta.rvs(size = 10000)
plt.hist(r, range=(a, b), density=True)
plt.show()

И графики

выборка