# 给定概率密度，生成随机数 python实现

17 篇文章 5 订阅

1. rejection sampling
2. invert the cdf
3. Metropolis Algorithm (MCMC)

PDF是概率分布函数，对其积分或者求和可以得到CDF（累积概率分布函数），PDF积分或求和的结果始终为1

1、根据pdf得到cdf
2、由cdf得到inverse of the cdf
3、对于给定的均匀分布[0,1),带入inverse cdf，得到的结果即是我们需要的x

1、当CDF的逆函数好求时，直接根据公式求取，
2、反之当CDF的逆函数不好求时，用数值模拟方法

import numpy as np
import math
import random
import matplotlib.pyplot as plt
import collections

count_dict = dict()
bin_count = 20

def inverseCDF():
"""
return the x value in PDF
"""
uniform_random = random.random()
return inverse_cdf(uniform_random)

def pdf(x):
return 2 * x

# cdf = x^2, 其逆函数很好求，因此直接用公式法
def inverse_cdf(x):
return math.sqrt(x)

def draw_pdf(D):
global bin_count
D = collections.OrderedDict(sorted(D.items()))
plt.bar(range(len(D)), list(D.values()), align='center')
# 因为映射bin的时候采用的floor操作，因此加上0.5
value_list = [(key + 0.5) / bin_count for key in D.keys()]
plt.xticks(range(len(D)), value_list)
plt.xlabel('x', fontsize=5)
plt.ylabel('counts', fontsize=5)
plt.title('counting bits')
plt.show()

for i in range(90000):
x = inverseCDF()
# 用bin去映射，否则不好操作
bin = math.floor(x * bin_count)  # type(bin): int
count_dict[bin] = count_dict.get(bin, 0) + 1

draw_pdf(count_dict)


table的size越大则结果越真实（即区间划分的个数）

import numpy as np
import math
import random
import matplotlib.pyplot as plt
import collections

lookup_table_size = 40
CDFlookup_table = np.zeros((lookup_table_size))

count_dict = dict()
bin_count = 20

def inverse_cdf_numerically(y):
global lookup_table_size
global CDFlookup_table
value = 0.0
for i in range(lookup_table_size):
x = i * 1.0 / (lookup_table_size - 1)
value += pdf2(x)
CDFlookup_table[i] = value
CDFlookup_table /= value  # normalize the cdf

if y < CDFlookup_table[0]:
t = y / CDFlookup_table[0]
return t / lookup_table_size
index = -1
for j in range(lookup_table_size):
if CDFlookup_table[j] >= y:
index = j
break
# linear interpolation
t = (y - CDFlookup_table[index - 1]) / \
(CDFlookup_table[index] - CDFlookup_table[index - 1])
fractional_index = index + t  # 因为index从0开始,所以不是 (index-1)+t
return fractional_index / lookup_table_size

def inverseCDF():
"""
return the x value in PDF
"""
uniform_random = random.random()
return inverse_cdf_numerically(uniform_random)

def pdf2(x):
return (x * x * x - 10.0 * x * x + 5.0 * x + 11.0) / (10.417)

def draw_pdf(D):
global bin_count
D = collections.OrderedDict(sorted(D.items()))
plt.bar(range(len(D)), list(D.values()), align='center')
value_list = [(key + 0.5) / bin_count for key in D.keys()]
plt.xticks(range(len(D)), value_list)
plt.xlabel('x', fontsize=5)
plt.ylabel('counts', fontsize=5)
plt.title('counting bits')
plt.show()

for i in range(90000):
x = inverseCDF()
bin = math.floor(x * bin_count)  # type(bin): int
count_dict[bin] = count_dict.get(bin, 0) + 1

draw_pdf(count_dict)


import numpy as np
import matplotlib.pyplot as plt
"""
reference:
https://blog.demofox.org/2017/07/25/counting-bits-the-normal-distribution/
"""

def plot_bar_x():
# this is for plotting purpose
index = np.arange(counting.shape[0])
plt.bar(index, counting)
plt.xlabel('x', fontsize=5)
plt.ylabel('counts', fontsize=5)
plt.title('counting bits')
plt.show()

# if dice_side=2, is binomial distribution
# if dice_side>2 , is multinomial distribution
dice_side = 2
# if N becomes larger, then multinomial distribution will more like normal distribution
N = 100

counting = np.zeros(((dice_side - 1) * N + 1))

for i in range(30000):
sum = 0
for j in range(N):
dice_result = np.random.randint(0, dice_side)
sum += dice_result

counting[sum] += 1

# normalization
counting /= np.sum(counting)
plot_bar_x()


• 11
点赞
• 69
收藏
觉得还不错? 一键收藏
• 3
评论
01-06 5489
03-24 2221
02-17 714
03-03 5218

### “相关推荐”对你有帮助么？

• 非常没帮助
• 没帮助
• 一般
• 有帮助
• 非常有帮助

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。