| 1 |
# svpmc: Stochastic Volatility Inference via Population Monte Carlo |
|---|
| 2 |
# |
|---|
| 3 |
# Copyright (C) 2007-2008 by Edwin A. Suominen, http://www.eepatents.com |
|---|
| 4 |
# |
|---|
| 5 |
# This program is free software; you can redistribute it and/or modify it under |
|---|
| 6 |
# the terms of the GNU General Public License as published by the Free Software |
|---|
| 7 |
# Foundation; either version 2 of the License, or (at your option) any later |
|---|
| 8 |
# version. |
|---|
| 9 |
# |
|---|
| 10 |
# This program is distributed in the hope that it will be useful, but WITHOUT |
|---|
| 11 |
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS |
|---|
| 12 |
# FOR A PARTICULAR PURPOSE. See the file COPYING for more details. |
|---|
| 13 |
# |
|---|
| 14 |
# You should have received a copy of the GNU General Public License along with |
|---|
| 15 |
# this program; if not, write to the Free Software Foundation, Inc., 51 |
|---|
| 16 |
# Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA |
|---|
| 17 |
|
|---|
| 18 |
""" |
|---|
| 19 |
Non-Uniform Sampling |
|---|
| 20 |
""" |
|---|
| 21 |
|
|---|
| 22 |
import scipy as s |
|---|
| 23 |
from scipy import random |
|---|
| 24 |
|
|---|
| 25 |
from weave import Weaver |
|---|
| 26 |
|
|---|
| 27 |
|
|---|
| 28 |
class Resampler(Weaver): |
|---|
| 29 |
""" |
|---|
| 30 |
Call an instance of me with a 1-D array of weights to get an index array |
|---|
| 31 |
that resamples values from some external array. Each element of the array |
|---|
| 32 |
will have a probability of being sampled (with replacement) that is |
|---|
| 33 |
proportional to its corresponding weight in the array. |
|---|
| 34 |
|
|---|
| 35 |
Uses Walker's alias algorithm as set forth in L. Devroye, Non-Uniform |
|---|
| 36 |
Random Variate Generation, p. 109, |
|---|
| 37 |
U{http://cg.scs.carleton.ca/~luc/rnbookindex.html}. |
|---|
| 38 |
""" |
|---|
| 39 |
|
|---|
| 40 |
def __init__(self, logWeights=False): |
|---|
| 41 |
self.logWeights = logWeights |
|---|
| 42 |
|
|---|
| 43 |
def __call__(self, W, N=None): |
|---|
| 44 |
""" |
|---|
| 45 |
If N is not specified, it defaults to the length of W, i.e., importance |
|---|
| 46 |
resampling. |
|---|
| 47 |
|
|---|
| 48 |
If my I{logWeights} attribute is set C{True}, the weights are |
|---|
| 49 |
transformed from log to linear. |
|---|
| 50 |
""" |
|---|
| 51 |
if self.logWeights: |
|---|
| 52 |
W = s.exp(W - W.max()) |
|---|
| 53 |
if not W.sum(): |
|---|
| 54 |
raise ValueError("Undefined result with all zero weights") |
|---|
| 55 |
if not s.isfinite(W).all(): |
|---|
| 56 |
W = 1 - s.isfinite(W) |
|---|
| 57 |
W = W.astype(float) / sum(W) |
|---|
| 58 |
I0 = (s.isfinite(W) * s.greater(W, 1E-7)).nonzero()[0] |
|---|
| 59 |
K = len(I0) |
|---|
| 60 |
if N is None: |
|---|
| 61 |
N = K |
|---|
| 62 |
if K == 0: |
|---|
| 63 |
return [] |
|---|
| 64 |
if K == 1: |
|---|
| 65 |
return I0.repeat(N) |
|---|
| 66 |
x = s.zeros((K,2)) |
|---|
| 67 |
x[:,0] = K*W[I0] |
|---|
| 68 |
a = s.less(x[:,0], 1.0).nonzero()[0] |
|---|
| 69 |
self.inline('x', 'a', 'b', b=s.setdiff1d(s.arange(0, K), a)) |
|---|
| 70 |
V = s.rand(N) |
|---|
| 71 |
RI = s.random.randint(0, K, N) |
|---|
| 72 |
I1 = s.greater(V, x[RI,0]).choose(RI, x[RI,1].astype(int)) |
|---|
| 73 |
return I0[I1] |
|---|
