Source Code for Module csb.statistics.ars

  1  """ 
  2  Adaptive Rejection Sampling (ARS) 
  3   
  4  The ARS class generates a single random sample from a 
  5  univariate distribution specified by an instance of the 
  6  LogProb class, implemented by the user. An instance of 
  7  LogProb returns the log of the probability density and 
  8  its derivative. The log probability function passed must 
  9  be concave. 
 10   
 11  The user must also supply initial guesses.  It is not 
 12  essential that these values be very accurate, but performance 
 13  will generally depend on their accuracy. 
 14  """ 
 15   
 16  from numpy import exp, log 
 17   
 18 -class Envelope(object): 
 19      """ 
 20      Envelope function for adaptive rejection sampling. 
 21       
 22      The envelope defines a piecewise linear upper and lower 
 23      bounding function of the concave log-probability. 
 24      """ 
 25 -    def __init__(self, x, h, dh): 
 26   
 27          from numpy import array, inf 
 28   
 29          self.x = array(x) 
 30          self.h = array(h) 
 31          self.dh = array(dh) 
 32          self.z0 = -inf 
 33          self.zk = inf 
 34           
 35 -    def z(self): 
 36          """ 
 37          Support intervals for upper bounding function. 
 38          """ 
 39          from numpy import concatenate 
 40   
 41          h = self.h 
 42          dh = self.dh 
 43          x = self.x 
 44   
 45          z = (h[1:] - h[:-1] + x[:-1] * dh[:-1] - x[1:] * dh[1:]) / \ 
 46              (dh[:-1] - dh[1:]) 
 47   
 48          return concatenate(([self.z0], z, [self.zk])) 
 49   
 50 -    def u(self, x): 
 51          """ 
 52          Piecewise linear upper bounding function. 
 53          """ 
 54          z = self.z()[1:-1] 
 55          j = (x > z).sum() 
 56   
 57          return self.h[j] + self.dh[j] * (x - self.x[j]) 
 58   
 59 -    def u_max(self): 
 60   
 61          z = self.z()[1:-1] 
 62   
 63          return (self.h + self.dh * (z - self.x)).max() 
 64   
 65 -    def l(self, x): 
 66          """ 
 67          Piecewise linear lower bounding function. 
 68          """ 
 69          from numpy import inf 
 70   
 71          j = (x > self.x).sum() 
 72   
 73          if j == 0 or j == len(self.x): 
 74              return -inf 
 75          else: 
 76              j -= 1 
 77              return ((self.x[j + 1] - x) * self.h[j] + (x - self.x[j]) * self.h[j + 1]) / \ 
 78                     (self.x[j + 1] - self.x[j]) 
 79   
 80 -    def insert(self, x, h, dh): 
 81          """ 
 82          Insert new support point for lower bounding function 
 83          (and indirectly for upper bounding function). 
 84          """ 
 85          from numpy import concatenate 
 86   
 87          j = (x > self.x).sum() 
 88   
 89          self.x = concatenate((self.x[:j], [x], self.x[j:])) 
 90          self.h = concatenate((self.h[:j], [h], self.h[j:])) 
 91          self.dh = concatenate((self.dh[:j], [dh], self.dh[j:])) 
 92   
 93 -    def log_masses(self): 
 94           
 95          from numpy import  abs, putmask 
 96   
 97          z = self.z() 
 98          b = self.h - self.x * self.dh 
 99          a = abs(self.dh) 
100          m = (self.dh > 0) 
101          q = self.x * 0.         
102          putmask(q, m, z[1:]) 
103          putmask(q, 1 - m, z[:-1]) 
104           
105          log_M = b - log(a) + log(1 - exp(-a * (z[1:] - z[:-1]))) + \ 
106                  self.dh * q 
107   
108          return log_M 
109   
110 -    def masses(self): 
111   
112          z = self.z() 
113          b = self.h - self.x * self.dh 
114          a = self.dh 
115           
116          return exp(b) * (exp(a * z[1:]) - exp(a * z[:-1])) / a 
117   
118 -    def sample(self): 
119   
120          from numpy.random import random 
121          from numpy import add 
122          from csb.numeric import log_sum_exp 
123           
124          log_m = self.log_masses() 
125          log_M = log_sum_exp(log_m) 
126          c = add.accumulate(exp(log_m - log_M)) 
127          u = random() 
128          j = (u > c).sum() 
129   
130          a = self.dh[j] 
131          z = self.z() 
132           
133          xmin, xmax = z[j], z[j + 1] 
134   
135          u = random() 
136   
137          if a > 0: 
138              return xmax + log(u + (1 - u) * exp(-a * (xmax - xmin))) / a 
139          else: 
140              return xmin + log(u + (1 - u) * exp(a * (xmax - xmin))) / a 
141   
142   
143 -class LogProb(object): 
144   
145 -    def __call__(self, x): 
146          raise NotImplementedError() 
147   
148 -class Gauss(LogProb): 
149   
150 -    def __init__(self, mu, sigma=1.): 
151   
152          self.mu = float(mu) 
153          self.sigma = float(sigma) 
154   
155 -    def __call__(self, x): 
156   
157          return -0.5 * (x - self.mu) ** 2 / self.sigma ** 2, \ 
158                 - (x - self.mu) / self.sigma ** 2 
159   
160   
161 -class ARS(object): 
162   
163      from numpy import inf 
164   
165 -    def __init__(self, logp): 
166   
167          self.logp = logp 
168   
169 -    def initialize(self, x, z0=-inf, zmax=inf): 
170   
171          from numpy import array 
172   
173          self.hull = Envelope(array(x), *self.logp(array(x))) 
174          self.hull.z0 = z0 
175          self.hull.zk = zmax 
176   
177 -    def sample(self, maxiter=100): 
178   
179          from numpy.random import random 
180   
181          for i in range(maxiter): 
182   
183              x = self.hull.sample() 
184              l = self.hull.l(x) 
185              u = self.hull.u(x) 
186              w = random() 
187   
188              if w <= exp(l - u): return x 
189   
190              h, dh = self.logp(x) 
191   
192              if w <= exp(h - u): return x 
193   
194              self.hull.insert(x, h, dh) 
195