样本和概率密度函数的转换

一、按概率密度函数生成样本

参照python – How to simulate from an (arbitrary) continuous probability distribution? – Stack Overflow

方法一（龟速）:

from scipy import stats
import numpy as np
class norm_distribution(stats.rv_continuous):
    def _pdf(self, x):
        return stats.norm.pdf(x, 0, 1) #输入pdf
distribution = norm_distribution(a=-2,b=2) #设置边界缩短运行时间
randdist1=distribution.rvs(size=1000)

方法二：

def randdist(x, pdf, nvals):
    """Produce nvals random samples from pdf(x), assuming constant spacing in x."""

    # get cumulative distribution from 0 to 1
    cumpdf = np.cumsum(pdf)
    cumpdf *= 1/cumpdf[-1]

    # input random values
    randv = np.random.uniform(size=nvals)

    # find where random values would go
    idx1 = np.searchsorted(cumpdf, randv)
    # get previous value, avoiding division by zero below
    idx0 = np.where(idx1==0, 0, idx1-1)
    idx1[idx0==0] = 1

    # do linear interpolation in x
    frac1 = (randv - cumpdf[idx0]) / (cumpdf[idx1] - cumpdf[idx0])
    randdist = x[idx0]*(1-frac1) + x[idx1]*frac1
    return randdist
x=np.linspace(-2,2,1000)
y=stats.norm.pdf(x, 0, 1)#输入pdf
randdist2=randdist(x, y, 1000)

二、已知样本, 求概率密度

接上面的代码,

from scipy.stats import gaussian_kde
import matplotlib.pyplot as plt
kde1 = gaussian_kde(randdist1)
y1 = kde1(x)
plt.plot(x, y1,label="stats.rv_continuous",lw=2)
kde2 = gaussian_kde(randdist2)
y2 = kde2(x)
plt.plot(x, y2,label="interpolation",lw=2)
randdist3=np.random.normal(loc=0.0, scale=1.0, size=1000)
kde3 = gaussian_kde(randdist3)
y3=kde3(x)
plt.plot(x, y3,label="random.normal",lw=2)
plt.plot(x, y,label="real",lw=2)
plt.legend()
plt.show()