我想创建一个基于 Levy 截断法则的自定义分布,它显示为
p(r) = (r + r0)**(-beta)*exp(-r/k)
。
所以我是这样定义的:
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)
假设我想找到 r = 0
和 r = 50km
之间的距离分布。然后:
nmin = 0
nmax = 50
my_cv = LevyPDF(a=nmin, b=nmax, name='LevyPDF')
x = np.linspace(nmin, nmax, (nmax-nmin)*2)
我不明白为什么:
sum(my_cv.cdf(x)) = 2.22
而不是 1
。
那么如何根据我定义的分布定义 N = 2000000
随机距离的直方图?
最佳答案
使用您的最小示例(稍微改编):
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')
要从您的随机变量中抽样,请使用 rvs()
方法来自 rv_continuous
类:
N = 50000
X = my_cv.rvs(size=N, random_state=1)
将返回一个大小为 (N,)
的数组,其中包含从您的分布中采样的随机变量。使用 random_state
选项来卡住您的示例并使您的脚本可重复(它为您的采样定义了随机种子)。
请注意,随着 N
缓慢增加,计算时间会急剧增加。
要绘制直方图,请使用 matplotlib
库,请参阅 hist
:
fig, axe = plt.subplots()
n, bins, patches = axe.hist(X, 50, normed=1, facecolor='green', alpha=0.75)
plt.show(axe)
下面是一个 40 自由度的卡方采样示例:
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)
它导致:
关于Python:如何定义自定义分布?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/46055690/