In [503]:

    

from __future__ import division

import math as m
import numpy as np
from scipy.stats import norm, uniform
from scipy import linalg

import pymstk.bridgesampling as bs
reload(bs)
    
class Model(object):
    def __init__(self, a=1):
        self.a   = a
        self.dist = norm(2,1)
        
    def rvs(self, size=1):
        return self.dist.rvs((size,1))
    
    def pdf(self, x):
        return self.a*self.dist.pdf(x)
    
    def logpdf(self, x):
        return log(self.a) + log(self.dist.pdf(x))


    

In [520]:

    

m = Model(20)
m_samples = m.rvs(2000)
m_lvalues = log(m.pdf(m_samples))

mvu = bs.MVU([[-5,5]])
mvn = bs.MVN([2],[1.5])


    

In [521]:

    

#%%timeit
z = bs.bs(m_samples, m_lvalues, m.logpdf, mvu, n2=400, niter=10, return_all=True, guess=log(1))
print exp(z)


    

    




[  1.          18.41866155  18.34100396  18.34159075  18.34158632
  18.34158635  18.34158635  18.34158635  18.34158635  18.34158635]

In [525]:

    

#%%timeit
z = bs.bs(m_samples, m_lvalues, m.logpdf, mvn, n2=400, niter=10, return_all=True, guess=log(1))
print exp(z)


    

    




[  1.          19.49967986  19.70322954  19.7044629   19.70447034
  19.70447039  19.70447039  19.70447039  19.70447039  19.70447039]

In [526]:

    

zu = [bs.bs(m_samples, m_lvalues, m.logpdf, mvu,  n2=1000, niter=10, return_all=False, guess=log(1)) for i in range(150)]
zn = [bs.bs(m_samples, m_lvalues, m.logpdf, mvn,  n2=1000, niter=10, return_all=False, guess=log(1)) for i in range(150)]


    

In [527]:

    

hist(exp(zu), fc='b', alpha=0.5, range=(18,22));
hist(exp(zn), fc='r', alpha=0.5, range=(18,22));


    

In [50]:

    

z1 = [bs.bs(m_samples, m_lvalues, m.logpdf, mvu, n2=2000, niter=10, return_all=False, guess=log(1)) for i in range(150)]
z2 = [bs.bs(m_samples, m_lvalues, m.logpdf, mvu, n2=4000, niter=10, return_all=False, guess=log(1)) for i in range(150)]


    

    




---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-50-8bf7cde81434> in <module>()
      1 z1 = [bs.bs(m_samples, m_lvalues, m.logpdf, mvu, n2=2000, niter=10, return_all=False, guess=log(1)) for i in range(150)]
----> 2 z2 = [bs.bs(m_samples, m_lvalues, m.logpdf, mvu, n2=4000, niter=10, return_all=False, guess=log(1)) for i in range(150)]

/home/parviainen/.local/lib/python2.7/site-packages/pymstk/bridgesampling.pyc in bs(x, q, q1fun, q2fun, means, sigmas, n2, niter, return_all, guess)
     36 
     37     lq11 = np.asarray(q)
---> 38     lq12 = np.array([q1fun(x2) for x2 in x2s])
     39     lq22 = np.array([q2fun(x2) for x2 in x2s])
     40     lq21 = np.array([q2fun(x1) for x1 in x1s])

<ipython-input-46-d59ad753dd2a> in logpdf(self, x)
     21 
     22     def logpdf(self, x):
---> 23         return log(self.a) + log(self.dist.pdf(x))
     24 
     25 

/usr/lib/python2.7/dist-packages/scipy/stats/distributions.pyc in pdf(self, x)
    333 
    334     def pdf(self, x):    #raises AttributeError in frozen discrete distribution
--> 335         return self.dist.pdf(x, *self.args, **self.kwds)
    336 
    337     def cdf(self, x):

/usr/lib/python2.7/dist-packages/scipy/stats/distributions.pyc in pdf(self, x, *args, **kwds)
   1108         output = zeros(shape(cond),'d')
   1109         putmask(output,(1-cond0)*array(cond1,bool),self.badvalue)
-> 1110         if any(cond):
   1111             goodargs = argsreduce(cond, *((x,)+args+(scale,)))
   1112             scale, goodargs = goodargs[-1], goodargs[:-1]

KeyboardInterrupt: 

In [ ]:

    

hist(exp(z1), fc='b', alpha=0.5, range=(1900,2200));
hist(exp(z2), fc='g', alpha=0.5, range=(1900,2200));


    

In [231]:

    

m = Model(2035.5)
m_samples = m.rvs(20000)
m_lvalues = m.logpdf(m_samples)

mvu = MVU([[-7,7]])


    

In [234]:

    

z1 = [bs.bs(m_samples, m_lvalues, m.logpdf, mvu, n2=2000, niter=10, return_all=False, guess=log(1)) for i in range(150)]
z2 = [bs.bs(m_samples, m_lvalues, m.logpdf, mvu, n2=10000, niter=10, return_all=False, guess=log(1)) for i in range(150)]


    

In [235]:

    

hist(exp(z1), fc='b', alpha=0.5, range=(1900,2200));
hist(exp(z2), fc='g', alpha=0.5, range=(1900,2200));