1# Copyright 2020 The PyMC Developers 2# 3# Licensed under the Apache License, Version 2.0 (the "License"); 4# you may not use this file except in compliance with the License. 5# You may obtain a copy of the License at 6# 7# http://www.apache.org/licenses/LICENSE-2.0 8# 9# Unless required by applicable law or agreed to in writing, software 10# distributed under the License is distributed on an "AS IS" BASIS, 11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 12# See the License for the specific language governing permissions and 13# limitations under the License. 14 15import numpy as np 16import numpy.random as nr 17import scipy.linalg 18import theano 19 20import pymc3 as pm 21 22from pymc3.distributions import draw_values 23from pymc3.step_methods.arraystep import ( 24 ArrayStep, 25 ArrayStepShared, 26 Competence, 27 PopulationArrayStepShared, 28 metrop_select, 29) 30from pymc3.theanof import floatX 31 32__all__ = [ 33 "Metropolis", 34 "DEMetropolis", 35 "DEMetropolisZ", 36 "BinaryMetropolis", 37 "BinaryGibbsMetropolis", 38 "CategoricalGibbsMetropolis", 39 "NormalProposal", 40 "CauchyProposal", 41 "LaplaceProposal", 42 "PoissonProposal", 43 "MultivariateNormalProposal", 44] 45 46# Available proposal distributions for Metropolis 47 48 49class Proposal: 50 def __init__(self, s): 51 self.s = s 52 53 54class NormalProposal(Proposal): 55 def __call__(self): 56 return nr.normal(scale=self.s) 57 58 59class UniformProposal(Proposal): 60 def __call__(self): 61 return nr.uniform(low=-self.s, high=self.s, size=len(self.s)) 62 63 64class CauchyProposal(Proposal): 65 def __call__(self): 66 return nr.standard_cauchy(size=np.size(self.s)) * self.s 67 68 69class LaplaceProposal(Proposal): 70 def __call__(self): 71 size = np.size(self.s) 72 return (nr.standard_exponential(size=size) - nr.standard_exponential(size=size)) * self.s 73 74 75class PoissonProposal(Proposal): 76 def __call__(self): 77 return nr.poisson(lam=self.s, size=np.size(self.s)) - self.s 78 79 80class MultivariateNormalProposal(Proposal): 81 def __init__(self, s): 82 n, m = s.shape 83 if n != m: 84 raise ValueError("Covariance matrix is not symmetric.") 85 self.n = n 86 self.chol = scipy.linalg.cholesky(s, lower=True) 87 88 def __call__(self, num_draws=None): 89 if num_draws is not None: 90 b = np.random.randn(self.n, num_draws) 91 return np.dot(self.chol, b).T 92 else: 93 b = np.random.randn(self.n) 94 return np.dot(self.chol, b) 95 96 97class Metropolis(ArrayStepShared): 98 """Metropolis-Hastings sampling step""" 99 100 name = "metropolis" 101 102 default_blocked = False 103 generates_stats = True 104 stats_dtypes = [ 105 { 106 "accept": np.float64, 107 "accepted": bool, 108 "tune": bool, 109 "scaling": np.float64, 110 } 111 ] 112 113 def __init__( 114 self, 115 vars=None, 116 S=None, 117 proposal_dist=None, 118 scaling=1.0, 119 tune=True, 120 tune_interval=100, 121 model=None, 122 mode=None, 123 **kwargs 124 ): 125 """Create an instance of a Metropolis stepper 126 127 Parameters 128 ---------- 129 vars: list 130 List of variables for sampler 131 S: standard deviation or covariance matrix 132 Some measure of variance to parameterize proposal distribution 133 proposal_dist: function 134 Function that returns zero-mean deviates when parameterized with 135 S (and n). Defaults to normal. 136 scaling: scalar or array 137 Initial scale factor for proposal. Defaults to 1. 138 tune: bool 139 Flag for tuning. Defaults to True. 140 tune_interval: int 141 The frequency of tuning. Defaults to 100 iterations. 142 model: PyMC Model 143 Optional model for sampling step. Defaults to None (taken from context). 144 mode: string or `Mode` instance. 145 compilation mode passed to Theano functions 146 """ 147 148 model = pm.modelcontext(model) 149 150 if vars is None: 151 vars = model.vars 152 vars = pm.inputvars(vars) 153 154 if S is None: 155 S = np.ones(sum(v.dsize for v in vars)) 156 157 if proposal_dist is not None: 158 self.proposal_dist = proposal_dist(S) 159 elif S.ndim == 1: 160 self.proposal_dist = NormalProposal(S) 161 elif S.ndim == 2: 162 self.proposal_dist = MultivariateNormalProposal(S) 163 else: 164 raise ValueError("Invalid rank for variance: %s" % S.ndim) 165 166 self.scaling = np.atleast_1d(scaling).astype("d") 167 self.tune = tune 168 self.tune_interval = tune_interval 169 self.steps_until_tune = tune_interval 170 self.accepted = 0 171 172 # Determine type of variables 173 self.discrete = np.concatenate( 174 [[v.dtype in pm.discrete_types] * (v.dsize or 1) for v in vars] 175 ) 176 self.any_discrete = self.discrete.any() 177 self.all_discrete = self.discrete.all() 178 179 # remember initial settings before tuning so they can be reset 180 self._untuned_settings = dict( 181 scaling=self.scaling, steps_until_tune=tune_interval, accepted=self.accepted 182 ) 183 184 self.mode = mode 185 186 shared = pm.make_shared_replacements(vars, model) 187 self.delta_logp = delta_logp(model.logpt, vars, shared) 188 super().__init__(vars, shared) 189 190 def reset_tuning(self): 191 """Resets the tuned sampler parameters to their initial values.""" 192 for attr, initial_value in self._untuned_settings.items(): 193 setattr(self, attr, initial_value) 194 return 195 196 def astep(self, q0): 197 if not self.steps_until_tune and self.tune: 198 # Tune scaling parameter 199 self.scaling = tune(self.scaling, self.accepted / float(self.tune_interval)) 200 # Reset counter 201 self.steps_until_tune = self.tune_interval 202 self.accepted = 0 203 204 delta = self.proposal_dist() * self.scaling 205 206 if self.any_discrete: 207 if self.all_discrete: 208 delta = np.round(delta, 0).astype("int64") 209 q0 = q0.astype("int64") 210 q = (q0 + delta).astype("int64") 211 else: 212 delta[self.discrete] = np.round(delta[self.discrete], 0) 213 q = q0 + delta 214 else: 215 q = floatX(q0 + delta) 216 217 accept = self.delta_logp(q, q0) 218 q_new, accepted = metrop_select(accept, q, q0) 219 self.accepted += accepted 220 221 self.steps_until_tune -= 1 222 223 stats = { 224 "tune": self.tune, 225 "scaling": self.scaling, 226 "accept": np.exp(accept), 227 "accepted": accepted, 228 } 229 230 return q_new, [stats] 231 232 @staticmethod 233 def competence(var, has_grad): 234 return Competence.COMPATIBLE 235 236 237def tune(scale, acc_rate): 238 """ 239 Tunes the scaling parameter for the proposal distribution 240 according to the acceptance rate over the last tune_interval: 241 242 Rate Variance adaptation 243 ---- ------------------- 244 <0.001 x 0.1 245 <0.05 x 0.5 246 <0.2 x 0.9 247 >0.5 x 1.1 248 >0.75 x 2 249 >0.95 x 10 250 251 """ 252 if acc_rate < 0.001: 253 # reduce by 90 percent 254 return scale * 0.1 255 elif acc_rate < 0.05: 256 # reduce by 50 percent 257 return scale * 0.5 258 elif acc_rate < 0.2: 259 # reduce by ten percent 260 return scale * 0.9 261 elif acc_rate > 0.95: 262 # increase by factor of ten 263 return scale * 10.0 264 elif acc_rate > 0.75: 265 # increase by double 266 return scale * 2.0 267 elif acc_rate > 0.5: 268 # increase by ten percent 269 return scale * 1.1 270 271 return scale 272 273 274class BinaryMetropolis(ArrayStep): 275 """Metropolis-Hastings optimized for binary variables 276 277 Parameters 278 ---------- 279 vars: list 280 List of variables for sampler 281 scaling: scalar or array 282 Initial scale factor for proposal. Defaults to 1. 283 tune: bool 284 Flag for tuning. Defaults to True. 285 tune_interval: int 286 The frequency of tuning. Defaults to 100 iterations. 287 model: PyMC Model 288 Optional model for sampling step. Defaults to None (taken from context). 289 290 """ 291 292 name = "binary_metropolis" 293 294 generates_stats = True 295 stats_dtypes = [ 296 { 297 "accept": np.float64, 298 "tune": bool, 299 "p_jump": np.float64, 300 } 301 ] 302 303 def __init__(self, vars, scaling=1.0, tune=True, tune_interval=100, model=None): 304 305 model = pm.modelcontext(model) 306 307 self.scaling = scaling 308 self.tune = tune 309 self.tune_interval = tune_interval 310 self.steps_until_tune = tune_interval 311 self.accepted = 0 312 313 if not all([v.dtype in pm.discrete_types for v in vars]): 314 raise ValueError("All variables must be Bernoulli for BinaryMetropolis") 315 316 super().__init__(vars, [model.fastlogp]) 317 318 def astep(self, q0, logp): 319 320 # Convert adaptive_scale_factor to a jump probability 321 p_jump = 1.0 - 0.5 ** self.scaling 322 323 rand_array = nr.random(q0.shape) 324 q = np.copy(q0) 325 # Locations where switches occur, according to p_jump 326 switch_locs = rand_array < p_jump 327 q[switch_locs] = True - q[switch_locs] 328 329 accept = logp(q) - logp(q0) 330 q_new, accepted = metrop_select(accept, q, q0) 331 self.accepted += accepted 332 333 stats = { 334 "tune": self.tune, 335 "accept": np.exp(accept), 336 "p_jump": p_jump, 337 } 338 339 return q_new, [stats] 340 341 @staticmethod 342 def competence(var): 343 """ 344 BinaryMetropolis is only suitable for binary (bool) 345 and Categorical variables with k=1. 346 """ 347 distribution = getattr(var.distribution, "parent_dist", var.distribution) 348 if isinstance(distribution, pm.Bernoulli) or (var.dtype in pm.bool_types): 349 return Competence.COMPATIBLE 350 elif isinstance(distribution, pm.Categorical) and (distribution.k == 2): 351 return Competence.COMPATIBLE 352 return Competence.INCOMPATIBLE 353 354 355class BinaryGibbsMetropolis(ArrayStep): 356 """A Metropolis-within-Gibbs step method optimized for binary variables 357 358 Parameters 359 ---------- 360 vars: list 361 List of variables for sampler 362 order: list or 'random' 363 List of integers indicating the Gibbs update order 364 e.g., [0, 2, 1, ...]. Default is random 365 transit_p: float 366 The diagonal of the transition kernel. A value > .5 gives anticorrelated proposals, 367 which resulting in more efficient antithetical sampling. Default is 0.8 368 model: PyMC Model 369 Optional model for sampling step. Defaults to None (taken from context). 370 371 """ 372 373 name = "binary_gibbs_metropolis" 374 375 def __init__(self, vars, order="random", transit_p=0.8, model=None): 376 377 model = pm.modelcontext(model) 378 379 # transition probabilities 380 self.transit_p = transit_p 381 382 self.dim = sum(v.dsize for v in vars) 383 384 if order == "random": 385 self.shuffle_dims = True 386 self.order = list(range(self.dim)) 387 else: 388 if sorted(order) != list(range(self.dim)): 389 raise ValueError("Argument 'order' has to be a permutation") 390 self.shuffle_dims = False 391 self.order = order 392 393 if not all([v.dtype in pm.discrete_types for v in vars]): 394 raise ValueError("All variables must be binary for BinaryGibbsMetropolis") 395 396 super().__init__(vars, [model.fastlogp]) 397 398 def astep(self, q0, logp): 399 order = self.order 400 if self.shuffle_dims: 401 nr.shuffle(order) 402 403 q = np.copy(q0) 404 logp_curr = logp(q) 405 406 for idx in order: 407 # No need to do metropolis update if the same value is proposed, 408 # as you will get the same value regardless of accepted or reject 409 if nr.rand() < self.transit_p: 410 curr_val, q[idx] = q[idx], True - q[idx] 411 logp_prop = logp(q) 412 q[idx], accepted = metrop_select(logp_prop - logp_curr, q[idx], curr_val) 413 if accepted: 414 logp_curr = logp_prop 415 416 return q 417 418 @staticmethod 419 def competence(var): 420 """ 421 BinaryMetropolis is only suitable for Bernoulli 422 and Categorical variables with k=2. 423 """ 424 distribution = getattr(var.distribution, "parent_dist", var.distribution) 425 if isinstance(distribution, pm.Bernoulli) or (var.dtype in pm.bool_types): 426 return Competence.IDEAL 427 elif isinstance(distribution, pm.Categorical) and (distribution.k == 2): 428 return Competence.IDEAL 429 return Competence.INCOMPATIBLE 430 431 432class CategoricalGibbsMetropolis(ArrayStep): 433 """A Metropolis-within-Gibbs step method optimized for categorical variables. 434 This step method works for Bernoulli variables as well, but it is not 435 optimized for them, like BinaryGibbsMetropolis is. Step method supports 436 two types of proposals: A uniform proposal and a proportional proposal, 437 which was introduced by Liu in his 1996 technical report 438 "Metropolized Gibbs Sampler: An Improvement". 439 """ 440 441 name = "categorical_gibbs_metropolis" 442 443 def __init__(self, vars, proposal="uniform", order="random", model=None): 444 445 model = pm.modelcontext(model) 446 vars = pm.inputvars(vars) 447 448 dimcats = [] 449 # The above variable is a list of pairs (aggregate dimension, number 450 # of categories). For example, if vars = [x, y] with x being a 2-D 451 # variable with M categories and y being a 3-D variable with N 452 # categories, we will have dimcats = [(0, M), (1, M), (2, N), (3, N), (4, N)]. 453 for v in vars: 454 distr = getattr(v.distribution, "parent_dist", v.distribution) 455 if isinstance(distr, pm.Categorical): 456 k = draw_values([distr.k])[0] 457 elif isinstance(distr, pm.Bernoulli) or (v.dtype in pm.bool_types): 458 k = 2 459 else: 460 raise ValueError( 461 "All variables must be categorical or binary" + "for CategoricalGibbsMetropolis" 462 ) 463 start = len(dimcats) 464 dimcats += [(dim, k) for dim in range(start, start + v.dsize)] 465 466 if order == "random": 467 self.shuffle_dims = True 468 self.dimcats = dimcats 469 else: 470 if sorted(order) != list(range(len(dimcats))): 471 raise ValueError("Argument 'order' has to be a permutation") 472 self.shuffle_dims = False 473 self.dimcats = [dimcats[j] for j in order] 474 475 if proposal == "uniform": 476 self.astep = self.astep_unif 477 elif proposal == "proportional": 478 # Use the optimized "Metropolized Gibbs Sampler" described in Liu96. 479 self.astep = self.astep_prop 480 else: 481 raise ValueError("Argument 'proposal' should either be 'uniform' or 'proportional'") 482 483 super().__init__(vars, [model.fastlogp]) 484 485 def astep_unif(self, q0, logp): 486 dimcats = self.dimcats 487 if self.shuffle_dims: 488 nr.shuffle(dimcats) 489 490 q = np.copy(q0) 491 logp_curr = logp(q) 492 493 for dim, k in dimcats: 494 curr_val, q[dim] = q[dim], sample_except(k, q[dim]) 495 logp_prop = logp(q) 496 q[dim], accepted = metrop_select(logp_prop - logp_curr, q[dim], curr_val) 497 if accepted: 498 logp_curr = logp_prop 499 return q 500 501 def astep_prop(self, q0, logp): 502 dimcats = self.dimcats 503 if self.shuffle_dims: 504 nr.shuffle(dimcats) 505 506 q = np.copy(q0) 507 logp_curr = logp(q) 508 509 for dim, k in dimcats: 510 logp_curr = self.metropolis_proportional(q, logp, logp_curr, dim, k) 511 512 return q 513 514 def metropolis_proportional(self, q, logp, logp_curr, dim, k): 515 given_cat = int(q[dim]) 516 log_probs = np.zeros(k) 517 log_probs[given_cat] = logp_curr 518 candidates = list(range(k)) 519 for candidate_cat in candidates: 520 if candidate_cat != given_cat: 521 q[dim] = candidate_cat 522 log_probs[candidate_cat] = logp(q) 523 probs = softmax(log_probs) 524 prob_curr, probs[given_cat] = probs[given_cat], 0.0 525 probs /= 1.0 - prob_curr 526 proposed_cat = nr.choice(candidates, p=probs) 527 accept_ratio = (1.0 - prob_curr) / (1.0 - probs[proposed_cat]) 528 if not np.isfinite(accept_ratio) or nr.uniform() >= accept_ratio: 529 q[dim] = given_cat 530 return logp_curr 531 q[dim] = proposed_cat 532 return log_probs[proposed_cat] 533 534 @staticmethod 535 def competence(var): 536 """ 537 CategoricalGibbsMetropolis is only suitable for Bernoulli and 538 Categorical variables. 539 """ 540 distribution = getattr(var.distribution, "parent_dist", var.distribution) 541 if isinstance(distribution, pm.Categorical): 542 if distribution.k > 2: 543 return Competence.IDEAL 544 return Competence.COMPATIBLE 545 elif isinstance(distribution, pm.Bernoulli) or (var.dtype in pm.bool_types): 546 return Competence.COMPATIBLE 547 return Competence.INCOMPATIBLE 548 549 550class DEMetropolis(PopulationArrayStepShared): 551 """ 552 Differential Evolution Metropolis sampling step. 553 554 Parameters 555 ---------- 556 lamb: float 557 Lambda parameter of the DE proposal mechanism. Defaults to 2.38 / sqrt(2 * ndim) 558 vars: list 559 List of variables for sampler 560 S: standard deviation or covariance matrix 561 Some measure of variance to parameterize proposal distribution 562 proposal_dist: function 563 Function that returns zero-mean deviates when parameterized with 564 S (and n). Defaults to Uniform(-S,+S). 565 scaling: scalar or array 566 Initial scale factor for epsilon. Defaults to 0.001 567 tune: str 568 Which hyperparameter to tune. Defaults to None, but can also be 'scaling' or 'lambda'. 569 tune_interval: int 570 The frequency of tuning. Defaults to 100 iterations. 571 model: PyMC Model 572 Optional model for sampling step. Defaults to None (taken from context). 573 mode: string or `Mode` instance. 574 compilation mode passed to Theano functions 575 576 References 577 ---------- 578 .. [Braak2006] Cajo C.F. ter Braak (2006). 579 A Markov Chain Monte Carlo version of the genetic algorithm 580 Differential Evolution: easy Bayesian computing for real parameter spaces. 581 Statistics and Computing 582 `link <https://doi.org/10.1007/s11222-006-8769-1>`__ 583 """ 584 585 name = "DEMetropolis" 586 587 default_blocked = True 588 generates_stats = True 589 stats_dtypes = [ 590 { 591 "accept": np.float64, 592 "accepted": bool, 593 "tune": bool, 594 "scaling": np.float64, 595 "lambda": np.float64, 596 } 597 ] 598 599 def __init__( 600 self, 601 vars=None, 602 S=None, 603 proposal_dist=None, 604 lamb=None, 605 scaling=0.001, 606 tune=None, 607 tune_interval=100, 608 model=None, 609 mode=None, 610 **kwargs 611 ): 612 613 model = pm.modelcontext(model) 614 615 if vars is None: 616 vars = model.cont_vars 617 vars = pm.inputvars(vars) 618 619 if S is None: 620 S = np.ones(model.ndim) 621 622 if proposal_dist is not None: 623 self.proposal_dist = proposal_dist(S) 624 else: 625 self.proposal_dist = UniformProposal(S) 626 627 self.scaling = np.atleast_1d(scaling).astype("d") 628 if lamb is None: 629 # default to the optimal lambda for normally distributed targets 630 lamb = 2.38 / np.sqrt(2 * model.ndim) 631 self.lamb = float(lamb) 632 if tune not in {None, "scaling", "lambda"}: 633 raise ValueError('The parameter "tune" must be one of {None, scaling, lambda}') 634 self.tune = tune 635 self.tune_interval = tune_interval 636 self.steps_until_tune = tune_interval 637 self.accepted = 0 638 639 self.mode = mode 640 641 shared = pm.make_shared_replacements(vars, model) 642 self.delta_logp = delta_logp(model.logpt, vars, shared) 643 super().__init__(vars, shared) 644 645 def astep(self, q0): 646 if not self.steps_until_tune and self.tune: 647 if self.tune == "scaling": 648 self.scaling = tune(self.scaling, self.accepted / float(self.tune_interval)) 649 elif self.tune == "lambda": 650 self.lamb = tune(self.lamb, self.accepted / float(self.tune_interval)) 651 # Reset counter 652 self.steps_until_tune = self.tune_interval 653 self.accepted = 0 654 655 epsilon = self.proposal_dist() * self.scaling 656 657 # differential evolution proposal 658 # select two other chains 659 ir1, ir2 = np.random.choice(self.other_chains, 2, replace=False) 660 r1 = self.bij.map(self.population[ir1]) 661 r2 = self.bij.map(self.population[ir2]) 662 # propose a jump 663 q = floatX(q0 + self.lamb * (r1 - r2) + epsilon) 664 665 accept = self.delta_logp(q, q0) 666 q_new, accepted = metrop_select(accept, q, q0) 667 self.accepted += accepted 668 669 self.steps_until_tune -= 1 670 671 stats = { 672 "tune": self.tune, 673 "scaling": self.scaling, 674 "lambda": self.lamb, 675 "accept": np.exp(accept), 676 "accepted": accepted, 677 } 678 679 return q_new, [stats] 680 681 @staticmethod 682 def competence(var, has_grad): 683 if var.dtype in pm.discrete_types: 684 return Competence.INCOMPATIBLE 685 return Competence.COMPATIBLE 686 687 688class DEMetropolisZ(ArrayStepShared): 689 """ 690 Adaptive Differential Evolution Metropolis sampling step that uses the past to inform jumps. 691 692 Parameters 693 ---------- 694 lamb: float 695 Lambda parameter of the DE proposal mechanism. Defaults to 2.38 / sqrt(2 * ndim) 696 vars: list 697 List of variables for sampler 698 S: standard deviation or covariance matrix 699 Some measure of variance to parameterize proposal distribution 700 proposal_dist: function 701 Function that returns zero-mean deviates when parameterized with 702 S (and n). Defaults to Uniform(-S,+S). 703 scaling: scalar or array 704 Initial scale factor for epsilon. Defaults to 0.001 705 tune: str 706 Which hyperparameter to tune. Defaults to 'lambda', but can also be 'scaling' or None. 707 tune_interval: int 708 The frequency of tuning. Defaults to 100 iterations. 709 tune_drop_fraction: float 710 Fraction of tuning steps that will be removed from the samplers history when the tuning ends. 711 Defaults to 0.9 - keeping the last 10% of tuning steps for good mixing while removing 90% of 712 potentially unconverged tuning positions. 713 model: PyMC Model 714 Optional model for sampling step. Defaults to None (taken from context). 715 mode: string or `Mode` instance. 716 compilation mode passed to Theano functions 717 718 References 719 ---------- 720 .. [Braak2006] Cajo C.F. ter Braak (2006). 721 Differential Evolution Markov Chain with snooker updater and fewer chains. 722 Statistics and Computing 723 `link <https://doi.org/10.1007/s11222-008-9104-9>`__ 724 """ 725 726 name = "DEMetropolisZ" 727 728 default_blocked = True 729 generates_stats = True 730 stats_dtypes = [ 731 { 732 "accept": np.float64, 733 "accepted": bool, 734 "tune": bool, 735 "scaling": np.float64, 736 "lambda": np.float64, 737 } 738 ] 739 740 def __init__( 741 self, 742 vars=None, 743 S=None, 744 proposal_dist=None, 745 lamb=None, 746 scaling=0.001, 747 tune="lambda", 748 tune_interval=100, 749 tune_drop_fraction: float = 0.9, 750 model=None, 751 mode=None, 752 **kwargs 753 ): 754 model = pm.modelcontext(model) 755 756 if vars is None: 757 vars = model.cont_vars 758 vars = pm.inputvars(vars) 759 760 if S is None: 761 S = np.ones(model.ndim) 762 763 if proposal_dist is not None: 764 self.proposal_dist = proposal_dist(S) 765 else: 766 self.proposal_dist = UniformProposal(S) 767 768 self.scaling = np.atleast_1d(scaling).astype("d") 769 if lamb is None: 770 # default to the optimal lambda for normally distributed targets 771 lamb = 2.38 / np.sqrt(2 * model.ndim) 772 self.lamb = float(lamb) 773 if tune not in {None, "scaling", "lambda"}: 774 raise ValueError('The parameter "tune" must be one of {None, scaling, lambda}') 775 self.tune = True 776 self.tune_target = tune 777 self.tune_interval = tune_interval 778 self.tune_drop_fraction = tune_drop_fraction 779 self.steps_until_tune = tune_interval 780 self.accepted = 0 781 782 # cache local history for the Z-proposals 783 self._history = [] 784 # remember initial settings before tuning so they can be reset 785 self._untuned_settings = dict( 786 scaling=self.scaling, 787 lamb=self.lamb, 788 steps_until_tune=tune_interval, 789 accepted=self.accepted, 790 ) 791 792 self.mode = mode 793 794 shared = pm.make_shared_replacements(vars, model) 795 self.delta_logp = delta_logp(model.logpt, vars, shared) 796 super().__init__(vars, shared) 797 798 def reset_tuning(self): 799 """Resets the tuned sampler parameters and history to their initial values.""" 800 # history can't be reset via the _untuned_settings dict because it's a list 801 self._history = [] 802 for attr, initial_value in self._untuned_settings.items(): 803 setattr(self, attr, initial_value) 804 return 805 806 def astep(self, q0): 807 # same tuning scheme as DEMetropolis 808 if not self.steps_until_tune and self.tune: 809 if self.tune_target == "scaling": 810 self.scaling = tune(self.scaling, self.accepted / float(self.tune_interval)) 811 elif self.tune_target == "lambda": 812 self.lamb = tune(self.lamb, self.accepted / float(self.tune_interval)) 813 # Reset counter 814 self.steps_until_tune = self.tune_interval 815 self.accepted = 0 816 817 epsilon = self.proposal_dist() * self.scaling 818 819 it = len(self._history) 820 # use the DE-MCMC-Z proposal scheme as soon as the history has 2 entries 821 if it > 1: 822 # differential evolution proposal 823 # select two other chains 824 iz1 = np.random.randint(it) 825 iz2 = np.random.randint(it) 826 while iz2 == iz1: 827 iz2 = np.random.randint(it) 828 829 z1 = self._history[iz1] 830 z2 = self._history[iz2] 831 # propose a jump 832 q = floatX(q0 + self.lamb * (z1 - z2) + epsilon) 833 else: 834 # propose just with noise in the first 2 iterations 835 q = floatX(q0 + epsilon) 836 837 accept = self.delta_logp(q, q0) 838 q_new, accepted = metrop_select(accept, q, q0) 839 self.accepted += accepted 840 self._history.append(q_new) 841 842 self.steps_until_tune -= 1 843 844 stats = { 845 "tune": self.tune, 846 "scaling": self.scaling, 847 "lambda": self.lamb, 848 "accept": np.exp(accept), 849 "accepted": accepted, 850 } 851 852 return q_new, [stats] 853 854 def stop_tuning(self): 855 """At the end of the tuning phase, this method removes the first x% of the history 856 so future proposals are not informed by unconverged tuning iterations. 857 """ 858 it = len(self._history) 859 n_drop = int(self.tune_drop_fraction * it) 860 self._history = self._history[n_drop:] 861 return super().stop_tuning() 862 863 @staticmethod 864 def competence(var, has_grad): 865 if var.dtype in pm.discrete_types: 866 return Competence.INCOMPATIBLE 867 return Competence.COMPATIBLE 868 869 870def sample_except(limit, excluded): 871 candidate = nr.choice(limit - 1) 872 if candidate >= excluded: 873 candidate += 1 874 return candidate 875 876 877def softmax(x): 878 e_x = np.exp(x - np.max(x)) 879 return e_x / np.sum(e_x, axis=0) 880 881 882def delta_logp(logp, vars, shared): 883 [logp0], inarray0 = pm.join_nonshared_inputs([logp], vars, shared) 884 885 tensor_type = inarray0.type 886 inarray1 = tensor_type("inarray1") 887 888 logp1 = pm.CallableTensor(logp0)(inarray1) 889 890 f = theano.function([inarray1, inarray0], logp1 - logp0) 891 f.trust_input = True 892 return f 893