slice.py 文件源码

def main(num_samples, burn, lag, w):
    alpha = 1.0
    N = 25
    Z = gu.simulate_crp(N, alpha)
    K = max(Z) + 1

    # CRP with gamma prior.
    log_pdf_fun = lambda alpha : gu.logp_crp_unorm(N, K, alpha) - alpha
    proposal_fun = lambda : np.random.gamma(1.0, 1.0)
    D = (0, float('Inf'))

    samples = su.slice_sample(proposal_fun, log_pdf_fun, D,
        num_samples=num_samples, burn=burn, lag=lag, w=w)

    minval = min(samples)
    maxval = max(samples)
    xvals = np.linspace(minval, maxval, 100)
    yvals = np.array([math.exp(log_pdf_fun(x)) for x in xvals])
    yvals /= trapz(xvals, yvals)

    fig, ax = plt.subplots(2,1)

    ax[0].hist(samples, 50, normed=True)

    ax[1].hist(samples, 100, normed=True)
    ax[1].plot(xvals,-yvals, c='red', lw=3, alpha=.8)
    ax[1].set_xlim(ax[0].get_xlim())
    ax[1].set_ylim(ax[0].get_ylim())
    plt.show()