如何将线性分数分布创建为自定义离散概率分布?
How to create a linear fractional distribution as a custom discrete probability distribution?
我定义了以下自定义概率分布:
import scipy.stats as st
# parameters
a = 3 / 16
b = 1
class linear_fractional(st.rv_discrete):
def _pdf(self, n):
if (n == 0):
return (a + b - 1) / (a + b)
else:
return (a * b ** (n - 1)) / (a + b) ** (n + 1)
LF = linear_fractional()
LF.rvs()
当我让脚本 运行 时,我收到一条冗长的错误消息:
Traceback (most recent call last):
File "C:/Users/thoma/PycharmProjects/Host_Parasite_Coevolution/Asymptotics.py", line 17, in <module> LF.rvs()
File "C:\Users\thoma\AppData\Local\Programs\Python\Python37-32\lib\site-packages\scipy\stats\_distn_infrastructure.py", line 2969, in rvs
return super(rv_discrete, self).rvs(*args, **kwargs)
...
RecursionError: maximum recursion depth exceeded while calling a Python object
如果我改为 LF.mean()
,我会得到
Fatal Python error: Cannot recover from stack overflow.
有谁知道这是为什么以及我该如何解决这个问题?我必须定义概率分布的上限吗?
根据 the docs and this post 中给出的示例,该方法需要进行一些修改。重要的是,由于它是离散分布,因此应使用 _pmf
而不是 _pdf
。此外,_pmf
将被 numpy 风格的数组调用 n
,n == 0
的工作方式完全不同。
当 n == 0
时 (a * b ** (n - 1)) / (a + b) ** (n + 1)
等于 (a + b - 1) / (a + b)
,我们可以对所有 n
使用第一个表达式。但是,当 b
是整数且 n = -1
时,numpy 会产生错误。将 b
与 1.0
相乘会将其更改为 numpy 不会给出此类错误的浮点数。如果多次使用相同的参数 a
和 b
,可能会生成冻结分布。
这是一个示例,它创建了生成样本的直方图,并将其与 pmf
进行比较。
import scipy.stats as st
import numpy as np
from matplotlib import pyplot as plt
class linear_fractional(st.rv_discrete):
def _pmf(self, n, a, b):
return (a * (1.0 * b) ** (n - 1)) / (a + b) ** (n + 1)
# parameters
a = 3 / 16
b = 1
LF = linear_fractional()
N = 10000
plt.hist(LF.rvs(a, b, size=N), bins=np.arange(-0.5, 50), ec='w', label='histogram of samples')
plt.plot(LF.pmf(np.arange(50), a, b) * N, 'ro', label='probability mass function (scaled)')
plt.legend(title=f'$a={a}; b={b}$')
plt.autoscale(enable=True, axis='x', tight=True)
plt.show()
LF.mean(a, b)
输出 5.33333333333286
散点图是说明分布样本的另一种方法:
plt.scatter(np.random.uniform(0, 1, N), LF.rvs(a, b, size=N), marker=',', alpha=0.2, lw=0, s=1, color='crimson')
PS:当b=1
时,此分布的公式等于geometric distribution加p = a/(a+1)
减1。这要快得多,因为它完全在 numpy 中计算。
samples = np.random.geometric(a/(a+1), size=1000) - 1
我定义了以下自定义概率分布:
import scipy.stats as st
# parameters
a = 3 / 16
b = 1
class linear_fractional(st.rv_discrete):
def _pdf(self, n):
if (n == 0):
return (a + b - 1) / (a + b)
else:
return (a * b ** (n - 1)) / (a + b) ** (n + 1)
LF = linear_fractional()
LF.rvs()
当我让脚本 运行 时,我收到一条冗长的错误消息:
Traceback (most recent call last):
File "C:/Users/thoma/PycharmProjects/Host_Parasite_Coevolution/Asymptotics.py", line 17, in <module> LF.rvs()
File "C:\Users\thoma\AppData\Local\Programs\Python\Python37-32\lib\site-packages\scipy\stats\_distn_infrastructure.py", line 2969, in rvs
return super(rv_discrete, self).rvs(*args, **kwargs)
...
RecursionError: maximum recursion depth exceeded while calling a Python object
如果我改为 LF.mean()
,我会得到
Fatal Python error: Cannot recover from stack overflow.
有谁知道这是为什么以及我该如何解决这个问题?我必须定义概率分布的上限吗?
根据 the docs and this post 中给出的示例,该方法需要进行一些修改。重要的是,由于它是离散分布,因此应使用 _pmf
而不是 _pdf
。此外,_pmf
将被 numpy 风格的数组调用 n
,n == 0
的工作方式完全不同。
当 n == 0
时 (a * b ** (n - 1)) / (a + b) ** (n + 1)
等于 (a + b - 1) / (a + b)
,我们可以对所有 n
使用第一个表达式。但是,当 b
是整数且 n = -1
时,numpy 会产生错误。将 b
与 1.0
相乘会将其更改为 numpy 不会给出此类错误的浮点数。如果多次使用相同的参数 a
和 b
,可能会生成冻结分布。
这是一个示例,它创建了生成样本的直方图,并将其与 pmf
进行比较。
import scipy.stats as st
import numpy as np
from matplotlib import pyplot as plt
class linear_fractional(st.rv_discrete):
def _pmf(self, n, a, b):
return (a * (1.0 * b) ** (n - 1)) / (a + b) ** (n + 1)
# parameters
a = 3 / 16
b = 1
LF = linear_fractional()
N = 10000
plt.hist(LF.rvs(a, b, size=N), bins=np.arange(-0.5, 50), ec='w', label='histogram of samples')
plt.plot(LF.pmf(np.arange(50), a, b) * N, 'ro', label='probability mass function (scaled)')
plt.legend(title=f'$a={a}; b={b}$')
plt.autoscale(enable=True, axis='x', tight=True)
plt.show()
LF.mean(a, b)
输出 5.33333333333286
散点图是说明分布样本的另一种方法:
plt.scatter(np.random.uniform(0, 1, N), LF.rvs(a, b, size=N), marker=',', alpha=0.2, lw=0, s=1, color='crimson')
PS:当b=1
时,此分布的公式等于geometric distribution加p = a/(a+1)
减1。这要快得多,因为它完全在 numpy 中计算。
samples = np.random.geometric(a/(a+1), size=1000) - 1