test_sampling.py 文件源码

def test_arms_via_lognormal (self):
        """ Test by sampling from a lognormal distribution."""

        N = 5000

        params = [ (0.0, 1.0, numpy.inf), (-2.0, 0.5, numpy.inf), (2.0, 2.0, 1000.0) ]

        for mu, sig, x_max in params:

            print "==========================\nSampling lognormal mean %.4f  sd %.4f" % (mu, sig)

            ln = distributions.LogNormal (mu, sig)
            s_ex = []
            while len(s_ex) < N:
                x = ln.sample (1)
                if x < x_max: s_ex.append (x)

            sampling.reset_statistics()

            s_ars = sampling.arms (ln.lpdf, ln.dx_lpdf, 0.0, x_max, n=N, num_points=5)

            pylab.hold(False)
            pylab.plot (s_ars)
            pylab.savefig('iter.png')

            pylab.plot (s_ars[1:N-1], s_ars[2:N], 'ro')
            pylab.savefig('lag1.png')

            s_ex.sort()
            s_ars.sort()

            stats = sampling.statistics()
            print "ARMS: Num MH rejections = %d" % (stats['mh_rejects'])

            self.assertEquals (N, stats['values_generated'])
            self.assertTrue (0 < stats['mh_rejects'])

            netutils.check_quantiles (self, s_ex, s_ars, N)