Python编程产生非均匀随机数的几种方法代码分享
作者:冬木远景 发布时间:2023-02-10 02:00:19
1.反变换法
设需产生分布函数为F(x)的连续随机数X。若已有[0,1]区间均匀分布随机数R,则产生X的反变换公式为:
F(x)=r, 即x=F-1(r)
反函数存在条件:如果函数y=f(x)是定义域D上的单调函数,那么f(x)一定有反函数存在,且反函数一定是单调的。分布函数F(x)为是一个单调递增函数,所以其反函数存在。从直观意义上理解,因为r一一对应着x,而在[0,1]均匀分布随机数R≤r的概率P(R≤r)=r。 因此,连续随机数X≤x的概率P(X≤x)=P(R≤r)=r=F(x)
即X的分布函数为F(x)。
例子:下面的代码使用反变换法在区间[0, 6]上生成随机数,其概率密度近似为P(x)=e-x
import numpy as np
import matplotlib.pyplot as plt
# probability distribution we're trying to calculate
p = lambda x: np.exp(-x)
# CDF of p
CDF = lambda x: 1-np.exp(-x)
# invert the CDF
invCDF = lambda x: -np.log(1-x)
# domain limits
xmin = 0 # the lower limit of our domain
xmax = 6 # the upper limit of our domain
# range limits
rmin = CDF(xmin)
rmax = CDF(xmax)
N = 10000 # the total of samples we wish to generate
# generate uniform samples in our range then invert the CDF
# to get samples of our target distribution
R = np.random.uniform(rmin, rmax, N)
X = invCDF(R)
# get the histogram info
hinfo = np.histogram(X,100)
# plot the histogram
plt.hist(X,bins=100, label=u'Samples');
# plot our (normalized) function
xvals=np.linspace(xmin, xmax, 1000)
plt.plot(xvals, hinfo[0][0]*p(xvals), 'r', label=u'p(x)')
# turn on the legend
plt.legend()
plt.show()
一般来说,直方图的外廓曲线接近于总体X的概率密度曲线。
2.舍选抽样法(Rejection Methold)
用反变换法生成随机数时,如果求不出F-1(x)的解析形式或者F(x)就没有解析形式,则可以用F-1(x)的近似公式代替。但是由于反函数计算量较大,有时也是很不适宜的。另一种方法是由Von Neumann提出的舍选抽样法。下图中曲线w(x)为概率密度函数,按该密度函数产生随机数的方法如下:
基本的rejection methold步骤如下:
1. Draw x uniformly from [xmin xmax]
2. Draw x uniformly from [0, ymax]
3. if y < w(x),accept the sample, otherwise reject it
4. repeat
即落在曲线w(x)和X轴所围成区域内的点接受,落在该区域外的点舍弃。
例子:下面的代码使用basic rejection sampling methold在区间[0, 10]上生成随机数,其概率密度近似为P(x)=e-x
# -*- coding: utf-8 -*-
'''
The following code produces samples that follow the distribution P(x)=e^−x
for x=[0, 10] and generates a histogram of the sampled distribution.
'''
import numpy as np
import matplotlib.pyplot as plt
P = lambda x: np.exp(-x)
# domain limits
xmin = 0 # the lower limit of our domain
xmax = 10 # the upper limit of our domain
# range limit (supremum) for y
ymax = 1
N = 10000 # the total of samples we wish to generate
accepted = 0 # the number of accepted samples
samples = np.zeros(N)
count = 0 # the total count of proposals
# generation loop
while (accepted < N):
# pick a uniform number on [xmin, xmax) (e.g. 0...10)
x = np.random.uniform(xmin, xmax)
# pick a uniform number on [0, ymax)
y = np.random.uniform(0,ymax)
# Do the accept/reject comparison
if y < P(x):
samples[accepted] = x
accepted += 1
count +=1
print count, accepted
# get the histogram info
# If bins is an int, it defines the number of equal-width bins in the given range
(n, bins)= np.histogram(samples, bins=30) # Returns: n-The values of the histogram,n是直方图中柱子的高度
# plot the histogram
plt.hist(samples,bins=30,label=u'Samples') # bins=30即直方图中有30根柱子
# plot our (normalized) function
xvals=np.linspace(xmin, xmax, 1000)
plt.plot(xvals, n[0]*P(xvals), 'r', label=u'P(x)')
# turn on the legend
plt.legend()
plt.show()
>>>
99552 10000
3.推广的舍取抽样法
从上图中可以看出,基本的rejection methold法抽样效率很低,因为随机数x和y是在区间[xmin xmax]和区间[0 ymax]上均匀分布的,产生的大部分点不会落在w(x)曲线之下(曲线e-x的形状一边高一边低,其曲线下的面积占矩形面积的比例很小,则舍选抽样效率很低)。为了改进简单舍选抽样法的效率,可以构造一个新的密度函数q(x)(called a proposal distribution from which we can readily draw samples),使它的形状接近p(x),并选择一个常数k使得kq(x)≥w(x)对于x定义域内的值都成立。对应下图,首先从分布q(z)中生成随机数z0,然后按均匀分布从区间[0 kq(z0)]生成一个随机数u0。 if u0 > p(z0) then the sample is rejected,otherwise u0 is retained. 即下图中灰色区域内的点都要舍弃。可见,由于随机点u0只出现在曲线kq(x)之下,且在q(x)较大处出现次数较多,从而大大提高了采样效率。显然q(x)形状越接近p(x),则采样效率越高。
根据上述思想,也可以表达采样规则如下:
1. Draw x from your proposal distribution q(x)
2. Draw y uniformly from [0 1]
3. if y < p(x)/kq(x) , accept the sample, otherwise reject it
4. repeat
下面例子中选择函数p(x)=1/(x+1)作为proposal distribution,k=1。曲线1/(x+1)的形状与e-x相近。
import numpy as np
import matplotlib.pyplot as plt
p = lambda x: np.exp(-x) # our distribution
g = lambda x: 1/(x+1) # our proposal pdf (we're choosing k to be 1)
CDFg = lambda x: np.log(x +1) # generates our proposal using inverse sampling
# domain limits
xmin = 0 # the lower limit of our domain
xmax = 10 # the upper limit of our domain
# range limits for inverse sampling
umin = CDFg(xmin)
umax = CDFg(xmax)
N = 10000 # the total of samples we wish to generate
accepted = 0 # the number of accepted samples
samples = np.zeros(N)
count = 0 # the total count of proposals
# generation loop
while (accepted < N):
# Sample from g using inverse sampling
u = np.random.uniform(umin, umax)
xproposal = np.exp(u) - 1
# pick a uniform number on [0, 1)
y = np.random.uniform(0, 1)
# Do the accept/reject comparison
if y < p(xproposal)/g(xproposal):
samples[accepted] = xproposal
accepted += 1
count +=1
print count, accepted
# get the histogram info
hinfo = np.histogram(samples,50)
# plot the histogram
plt.hist(samples,bins=50, label=u'Samples');
# plot our (normalized) function
xvals=np.linspace(xmin, xmax, 1000)
plt.plot(xvals, hinfo[0][0]*p(xvals), 'r', label=u'p(x)')
# turn on the legend
plt.legend()
plt.show()
>>>
24051 10000
可以对比基本的舍取法和改进的舍取法的结果,前者产生符合要求分布的10000个随机数运算了99552步,后者运算了24051步,可以看到效率明显提高。
总结
Python数据可视化编程通过Matplotlib创建散点图代码示例
Python编程实现使用线性回归预测数据
Python数据可视化正态分布简单分析及实现代码
如有不足之处,欢迎留言指出。感谢朋友们对本站的支持!
来源:http://www.cnblogs.com/21207-iHome/p/5266399.html


猜你喜欢
- 在使用柱状图时,经常遇到需要多组数据进行比较的情况。绘制单个数据系列的柱形图比较简单,多组数据柱状图绘制的关键有三点:多次调用bar()函数
- 最近好多伙伴说,我用vue做的项目本地是可以的,但部署到服务器遇到好多问题:资源找不到,直接访问index.html页面空白,刷新当前路由4
- 1.直接输入sql执行MySQL> select now();+---------------------+| now() |+---
- 本文实例讲述了Python实现根据日期获取当天凌晨时间戳的方法。分享给大家供大家参考,具体如下:# -*- coding:utf-8 -*-
- 如下所示:#coding:utf-8 ''''' Created on 2014-7-24 @aut
- 本文记录了mysql 5.7.16安装配置方法,具体内容如下第一步:下载下载地址滚动到下方就能看到了,根据自己的需求下载;我的电脑为64为的
- 引用body标签有两做法: 第一种:使用DOM Core 即引用某个给定文档的第一个(也是仅有的一个)body标签 document.get
- 目录结构:只需在自己的python项目下随便创建一个文件夹(下图中为:daka),然后将下载的chromedriver.exe、ask_fo
- qqbot 是一个用 python 实现的、基于腾讯 SmartQQ 协议的 QQ 机器人框架,可运行在 Linux 、 Windows 和
- 由于我们分发的python应用可能运行在64位环境,也可能运行在32位环境,所以我们需要为同一套应用代码配置两套打包环境,怎么配置?步骤如下
- 一、ZeroMQ概述 ZeroMQ(又名ØMQ,MQ,或zmq)像一个可嵌入的网络库,但其作用就像一个并发框
- 默认级别:warningimport logginglogging.debug('debug message')loggin
- Beautiful Soup使用时,一般可以通过指定对应的name和attrs去搜索,特定的名字和属性,以找到所需要的部分的html代码。但
- python opencv把一张图片嵌入(叠加)到另一张图片上1、背景:最近做了个烟火生成系统的界面设计,需要将烟雾图片嵌入到任意一张图片中
- Django中获取text,password名字:<input type="text" name="na
- 访问指定节点: getElementsByName(): <html> <head> <title>DO
- 1、查询语句的执行顺序select[distinct] from join(如left 
- SQL查询中什么时候需要使用表别名?今天写MySQL时遇到使用表别名的问题,这里重新总结一下。1、 表名很长时select * from w
- 任何一个行业里,当有一头近乎垄断的大象盘踞着的时候,生活在大象身后的蚂蚁们既是悲哀又是幸运的。悲哀的是市场已近乎被大象垄断留给他们的空间已经
- 相信为数不少的系统管理员每天都在做着同一样的工作——对数据进行备份。一旦哪一天疏忽了,而这一天系统又恰恰发生了故障,需要进行数据恢复,那么此