python - 如何生成具有预定义概率分布的随机数?
问题描述
我想在 python 中实现一个函数(使用numpy
),该函数将数学函数(例如,p(x) = e^(-x)
如下所示)作为输入并生成随机数,这些随机数根据该数学函数的概率分布进行分布。我需要绘制它们,这样我们才能看到分布。
我实际上需要一个随机数生成器函数,用于以下两个数学函数作为输入,但如果它可以采用其他函数,为什么不呢:
1) p(x) = e^(-x)
2)g(x) = (1/sqrt(2*pi)) * e^(-(x^2)/2)
有谁知道这在 python 中是如何实现的?
解决方案
对于您需要的简单分布,或者如果您有一个易于反转的封闭形式 CDF,您可以在 NumPy 中找到大量采样器,正如 Olivier 的回答中正确指出的那样。
对于任意分布,您可以使用马尔可夫链蒙特卡罗采样方法。
这些算法的最简单且可能更容易理解的变体是Metropolis采样。
基本思想是这样的:
- 从一个随机点开始,随机
x
走一步xnew = x + delta
- 评估起点
p(x)
和新起点的期望概率分布p(xnew)
- 如果新点更有可能
p(xnew)/p(x) >= 1
接受移动 - 如果新点不太可能随机决定是接受还是拒绝取决于新点的可能性1
- 从这一点开始新的步骤并重复循环
可以看到,例如Sokal 2,用这种方法采样的点遵循接受概率分布。
可以在PyMC3
包中找到 Python 中 Montecarlo 方法的广泛实现。
示例实现
这是一个玩具示例,只是为了向您展示基本思想,而不是以任何方式作为参考实现。任何严肃的工作请参考成熟的包。
def uniform_proposal(x, delta=2.0):
return np.random.uniform(x - delta, x + delta)
def metropolis_sampler(p, nsamples, proposal=uniform_proposal):
x = 1 # start somewhere
for i in range(nsamples):
trial = proposal(x) # random neighbour from the proposal distribution
acceptance = p(trial)/p(x)
# accept the move conditionally
if np.random.uniform() < acceptance:
x = trial
yield x
让我们看看它是否适用于一些简单的发行版
高斯混合
def gaussian(x, mu, sigma):
return 1./sigma/np.sqrt(2*np.pi)*np.exp(-((x-mu)**2)/2./sigma/sigma)
p = lambda x: gaussian(x, 1, 0.3) + gaussian(x, -1, 0.1) + gaussian(x, 3, 0.2)
samples = list(metropolis_sampler(p, 100000))
柯西
def cauchy(x, mu, gamma):
return 1./(np.pi*gamma*(1.+((x-mu)/gamma)**2))
p = lambda x: cauchy(x, -2, 0.5)
samples = list(metropolis_sampler(p, 100000))
任意函数
您实际上不必从适当的概率分布中进行抽样。您可能只需要强制执行一个有限域来对您的随机步骤进行采样3
p = lambda x: np.sqrt(x)
samples = list(metropolis_sampler(p, 100000, domain=(0, 10)))
p = lambda x: (np.sin(x)/x)**2
samples = list(metropolis_sampler(p, 100000, domain=(-4*np.pi, 4*np.pi)))
结论
关于提案分布、收敛性、相关性、效率、应用程序、贝叶斯形式主义、其他 MCMC 采样器等,还有很多话要说。我认为这不是合适的地方,还有很多比什么更好的材料我可以在这里写在线可用。
这里的想法是支持概率较高的探索,但仍然关注低概率区域,因为它们可能会导致其他峰值。基础是提案分布的选择,即你如何选择新的点来探索。太小的步骤可能会将您限制在分布的有限区域,太大可能会导致非常低效的探索。
面向物理。如今,贝叶斯形式主义(Metropolis-Hastings)是首选,但恕我直言,初学者更难掌握。网上有很多可用的教程,例如杜克大学的这个。
未显示实现不会增加太多混乱,但很简单,您只需在域边缘包装试验步骤或使所需功能在域外归零。
推荐阅读
- video-streaming - 重新定位的 Moov atom 不允许 chrome 播放视频,但 firefox 可以播放
- android - 在 sharedPreferences 中保存单个对象属性
- angular - 当插值连接“线性梯度”作为“背景图像”的值的条件时,图像消失。为什么?
- javascript - 在反应应用程序中包含调度程序模块
- windows - 将 Eclipse 复制到另一个用户的帐户时出现问题
- angular - Angular primeng 如何将带有主机的组件样式移动到通用样式
- remote-desktop - 如何设置任务计划程序任务以远程桌面连接到另一台机器?
- apache-nifi - 如何修复无法索引出处事件错误?
- javascript - 将 Javascript 变量数据粘贴到 SharePoint Online 列表表单中
- facebook - Facebook Ads 中应用延迟深度链接的生效时间是多少?