代码之家  ›  专栏  ›  技术社区  ›  user8188120

emcee MCMC采样未收敛到适当的参数值

  •  0
  • user8188120  · 技术社区  · 6 年前

    emcee 在Python中用一个预定义的似然函数进行MCMC采样,以找到两个数据总体之间的最佳边界。

    为了 http://dfm.io/emcee/current/user/line/

    在给定一些线性边界线的情况下,似然函数计算真正分类和真负分类,并用于最小化两个值之间的差异,同时最大化它们的总和。

    m b ,梯度和偏移(或偏差),对于边界线,我得到了一些非常大和/或非常小的值。

    我在下面放了一个示例代码,它生成一些划分得很好的总体,然后围绕参数值的初始猜测生成mcmc。我不确定为什么MCMC链不能很好地收敛到适当的值,所以任何帮助都将不胜感激。

    下面的代码应该是现成的。

    import emcee
    import numpy as np 
    from sklearn.metrics import confusion_matrix
    import matplotlib.pyplot as plt
    
    #generate some test x and y data
    folded_xy_train = np.random.uniform(0,1,10000) #test x data
    folded_z_train = np.random.uniform(0,1,10000) #test y data
    #define the true gradient and offset for the boundary line
    m_true, b_true = 5,-2.5
    
    #generate labels for the test data
    rounded_labels_train = np.ones(len(folded_z_train))
    model = (m_true*folded_xy_train) + b_true
    difference = model - folded_z_train
    rounded_labels_train[difference<0] = 0
    
    #show the test data
    plt.figure()
    plt.scatter(folded_xy_train,folded_z_train,c=rounded_labels_train,s=1.0)
    
    #define a likelihood function for the boundary line       
    def lnlike(theta, x, y, labels):
            m, b = theta
            model = (m*x) + b
            difference = model - y
            classifications = np.ones(len(y))
            classifications[difference<0]=0
            cfm = confusion_matrix(labels,classifications)
            cm = cfm.astype('float') / cfm.sum(axis=1)[:, np.newaxis]
            tn, fp, fn, tp = cm.ravel()
            likelihood_val = (0.5*(tp+tn))/(1+np.abs(tp-tn))
            ln_like = -np.log(likelihood_val)
            return ln_like
    
    #define a wide flat prior         
    def lnprior(theta):
        m, b, = theta
        if 0 < m < 10 and -20 < b < 5:
            return 0.0
        return -np.inf 
    
    #define the posterior         
    def lnprob(p, x, y, labels):
        lp = lnprior(p)
        if not np.isfinite(lp):
            return 0
        return lp + lnlike(p, x, y, labels)
    
    #setup the MCMC sampling  
    nwalkers = 4
    ndim = 2
    p0 = np.array([4.2,-2]) + [np.random.rand(ndim) for i in range(nwalkers)]
    sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(folded_xy_train, folded_z_train, rounded_labels_train))
    sampler.run_mcmc(p0, 500)
    
    #extract the MCMC paramater value chains       
    samples = sampler.chain[:, 50:, :].reshape((-1, ndim))        
    
    #view the parameter chains        
    plt.figure()
    plt.subplot(211)
    plt.plot(samples[:,0])
    plt.subplot(212)
    plt.plot(samples[:,1])     
    

    初始测试数据,显示给定x y数据的明显边界线(由二进制类标签着色):

    enter image description here

    样本移动,显示渐变参数(顶部)和偏移参数(底部)的奇怪采样。x轴表示MCMC行走步数,y轴表示给定步数下的MCMC参数值:

    enter image description here

    0 回复  |  直到 6 年前
    推荐文章