python实现gan代码_随机采样方法代码实现python

微信公众号：数据皮皮侠如果你觉得该公众号对你有帮助，欢迎关注、推广和宣传

内容目录：随机采样方法代码实现-python

文章目录1. 均匀分布，使用线性同余器2. Box-Muller变换正态分布3. 接受-拒绝采样4. MCMC采样变体-Metropolis-Hastings算法5. Gibbs抽样

文章目录

1. 均匀分布，使用线性同余器

from datetime import datetime
from math import sqrt, cos, pi, log
import matplotlib.pyplot as plt
import numpy as np
def uniform_gen(seed, a=0, b=1):
    _A = 214013
    _C = 214013
    _M = 2 ** 32
    while True:
        seed = (_A * seed + _C) % _M
        res = a + (b - a) * seed / _M
        yield res
def random_uniform(num, seed=None):
    res = []
    if not seed:
        seed = datetime.now().second
    gen = uniform_gen(seed)
    if num == 1:
        return gen.__next__()
    for i in range(num):
        res.append(gen.__next__())
    return np.array(res)
def plot_hist(x, bin=150, path = None):
    plt.hist(x, bin)
    if path:
        plt.savefig(path)
    plt.show()
res = random_uniform(10000)
plot_hist(res, 50, 'uniform.jpg')

2. Box-Muller变换正态分布

def box_muller(num, loc=0., scale=1.):
    seed1 = datetime.now().second
    seed2 = datetime.now().minute
    gen1 = uniform_gen(seed1)
    gen2 = uniform_gen(seed2)
    res = []
    for i in range(num):
        u1 = gen1.__next__()
        u2 = gen2.__next__()
        u = sqrt(-2 * log(u1)) * cos(2 * pi * u2)
        u = (u + loc) * scale
        res.append(u)
    return np.array(res)
res = box_muller(1000)
plot_hist(res, path='box_nuller.jpg')

3. 接受-拒绝采样

目标是生成服从代码块中的p(x)密度函数分布的样本，用kq(x)密度函数覆盖住p(x)，使得p(x)<=kq(x)。具体原理参考随机采样方法整理与讲解(MCMC、Gibbs Sampling等)不好的地方是很难找到合适的q分布，k的确定也麻烦。

def p(x):
    return (0.3 * np.exp(-(x - 0.3) ** 2) + 0.7 * np.exp(-(x - 2.) ** 2 / 0.3)) / 1.2113

def q(x, loc=0, scale=1):
    return 1 / (np.sqrt(2 * np.pi) * scale) * np.exp(-0.5 * (x - loc) ** 2 / scale ** 2)
def reject_sampling(pdf, k, loc, scale, num=1000):
    normal_sampling = box_muller(num, loc=loc, scale=scale)
    origin_pdf = pdf(normal_sampling)
    cover_pdf = q(normal_sampling, loc, scale) * k
    u = random_uniform(num) * cover_pdf
    res = normal_sampling[origin_pdf > u]
    return res
res = reject_sampling(p, 10, 1.4, 1.2, int(1e+7))
plot_hist(res)

4. MCMC采样变体-Metropolis-Hastings算法

def p(x):
    return 1 / (2 * np.pi) * np.exp(-x[0] ** 2 - x[1] ** 2)
def uniform_gen():
    # 计算定义域一个随机值
    return np.random.random() * 4 - 2
def metropolis(x0):
    x = (uniform_gen(), uniform_gen())
    u = np.random.uniform()
    alpha = min(1, p(x) / p(x0))
    if u         return x
    else:
        return x0
def gen_px(num):
    x = (uniform_gen(), uniform_gen())
    res = [x]
    for i in range(num):
        x = metropolis(x)
        res.append(x)
    return np.array(res)
def plt_scatter(x, path=None):
    x = np.array(x)
    plt.scatter(x[:, 0], x[:, 1])
    if path:
        plt.savefig(path)
    plt.show()
res = gen_px(5000)
plt_scatter(res, 'MCMC.jpg')

5. Gibbs抽样

def one_dim_transit(x, dim):
    candidates = []
    probs = []
    for i in range(10):
        tmp = x[:]
        tmp[dim] = uniform_gen()
        candidates.append(tmp[dim])
        probs.append(p(tmp))
    norm_probs = np.array(probs) / sum(probs)
    u = np.random.uniform()
    sum_probs = 0
    for i in range(10):
        sum_probs += norm_probs[i]
        if sum_probs >= u:
            return candidates[i]
def gibbs(x):
    dim = len(x)
    new = x[:]
    for i in range(dim):
        new[i] = one_dim_transit(new, i)
    return new
def gibbs_sampling(num=1000):
    x = [uniform_gen(), uniform_gen()]
    res = [x]
    for i in range(num):
        res.append(gibbs(x))
    return res    
res = gibbs_sampling(2000)
plt_scatter(res, path='gibbs.jpg')