Permalink
Browse files

refactoring. but something is not right

  • Loading branch information...
1 parent 544517b commit e68ac66e223caec9b16fe7edc21ac6ed3cdaf87e @wesm committed Mar 31, 2011
Showing with 410 additions and 108 deletions.
  1. +1 −1 examples/pw_ch07.py
  2. +13 −15 examples/wh_ch02.py
  3. +8 −7 examples/wh_ch03.py
  4. +6 −6 examples/wh_ch04.py
  5. +6 −4 examples/wh_ch08.py
  6. +6 −7 examples/wh_ch12.py
  7. +108 −9 statlib/arma.py
  8. +10 −1 statlib/datasets.py
  9. +245 −48 statlib/dlm.py
  10. +7 −10 statlib/multi.py
View
2 examples/pw_ch07.py
@@ -39,4 +39,4 @@ def simulate_76(pi=0.9, phi=0.9, v=1, T=5000):
yenusd = np.loadtxt('statlib/data/japan-usa1000.txt', delimiter=',')
ukusd = np.loadtxt('statlib/data/uk-usa1000.txt', delimiter=',')
-b
+
View
28 examples/wh_ch02.py
@@ -12,10 +12,10 @@
def ex_21():
y = datasets.table_22()
- x = np.ones(len(y), dtype=float)
+ x = [1]
- mean_prior = (0, 1)
- var_prior = (1, 0.01)
+ m0, C0 = (0, 1)
+ n0, s0 = (1, 0.01)
discount_factors = np.arange(0.6, 1.01, 0.05)
@@ -38,8 +38,8 @@ def ex_21():
mu = []
for disc in discount_factors:
- model = dlm.DLM(y, x, mean_prior=mean_prior,
- var_prior=var_prior, discount=disc)
+ model = dlm.DLM(y, x, m0=m0, C0=C0, n0=n0, s0=s0,
+ state_discount=disc)
ax1.plot(rng, model.forecast, '%.2f' % (1 - disc))
@@ -110,14 +110,12 @@ def boot_ci(sample, samples=10000):
y = datasets.table_22()
x = [1]
- mean_prior = (0, 1)
- var_prior = (1, 0.01)
-
- model3 = DLM(y, x, mean_prior=mean_prior, var_prior=var_prior,
- discount=0.1)
-
- model = DLM(y, x, mean_prior=mean_prior, var_prior=var_prior,
- discount=.9)
- model2 = DLM(y, x, mean_prior=mean_prior, var_prior=var_prior,
- discount=0.9)
+ m0, C0 = (0, 1)
+ n0, s0 = (1, 0.01)
+ model = DLM(y, x, m0=m0, C0=C0, n0=n0, s0=s0,
+ state_discount=.9)
+ model2 = DLM(y, x, m0=m0, C0=C0, n0=n0, s0=s0,
+ state_discount=0.9)
+ model3 = DLM(y, x, m0=m0, C0=C0, n0=n0, s0=s0,
+ state_discount=0.1)
View
15 examples/wh_ch03.py
@@ -22,15 +22,16 @@ def ex_310():
fig, axes = plt.subplots(3, 3, sharey=True, figsize=(12, 12))
for i, disc in enumerate(discount_factors):
- model = dlm.DLM(y, x, mean_prior=mean_prior,
- var_prior=var_prior, discount=disc)
+ model = dlm.DLM(y, x, m0=m0, C0=C0, n0=n0, s0=s0,
+ state_discount=disc)
rmse.append(model.rmse)
mad.append(model.mad)
like.append(model.pred_like)
var_est.append(model.var_est[-1])
- pred_var.append(model.forc_var[-1]) # model.var_est[-1] + model.mu_scale[-1] / disc
+ # model.var_est[-1] + model.mu_scale[-1] / disc
+ pred_var.append(model.forc_var[-1])
ax = axes[i / rows][i % rows]
@@ -62,8 +63,8 @@ def ex_310():
if __name__ == '__main__':
y, x = datasets.table_33()
- mean_prior = (0.45, 0.0025)
- var_prior = (1, 1)
+ m0, C0 = (0.45, 0.0025)
+ n0, s0 = (1, 1)
- model = dlm.DLM(y, x, mean_prior=mean_prior,
- var_prior=var_prior, discount=0.6)
+ model = dlm.DLM(y, x, m0=m0, C0=C0, n0=n0, s0=s0,
+ state_discount=0.6)
View
12 examples/wh_ch04.py
@@ -13,18 +13,18 @@
X = np.vstack((np.ones(len(x)),x )).T
-mean_prior = ([0, 0.45], [[0.005, 0],
+m0, C0 = ([0, 0.45], [[0.005, 0],
[0, 0.0025]])
-var_prior = (1, 1)
+n0, s0 = (1, 1)
G = np.diag([lam, phis[0]])
like = []
for phi in phis:
G = np.diag([lam, phi])
- model = DLM(y, X, G=G, mean_prior=mean_prior,
- var_prior=var_prior, discount=disc)
+ model = DLM(y, X, G=G, m0=m0, C0=C0, n0=n0, s0=s0,
+ state_discount=disc)
like.append(model.pred_like)
@@ -39,5 +39,5 @@
ax2.plot(phis, llr, label='LogLR')
plt.show()
-model = dlm.DLM(y, X, G=G, mean_prior=mean_prior,
- var_prior=var_prior, discount=disc)
+model = dlm.DLM(y, X, G=G, m0=m0, C0=C0, n0=n0, s0=s0,
+ state_discount=disc)
View
10 examples/wh_ch08.py
@@ -10,10 +10,12 @@
gas = datasets.table_81()
p = 12
-var_prior = (1, 0.05)
model = Polynomial(2) + FullEffectsFourier(12, harmonics=[1, 4])
# model = Polynomial(1) + FullEffectsFourier(12)
k = model.F.shape[1]
-mean_prior = (np.zeros(k), np.eye(k))
-dlm = DLM(gas, model.F, G=model.G, mean_prior=mean_prior,
- var_prior=var_prior, discount=.95)
+m0, C0 = (np.zeros(k), np.eye(k))
+n0, s0 = (1, 0.05)
+dlm = DLM(gas, model.F, G=model.G,
+ m0=m0, C0=C0, n0=n0, s0=s0,
+ state_discount=.95,
+ var_discount=1)
View
13 examples/wh_ch12.py
@@ -12,6 +12,7 @@
from statlib.dlm import *
from statlib.multi import *
import statlib.datasets as datasets
+import statlib.tools as tools
from pandas import DataMatrix
@@ -46,8 +47,7 @@ def get_multi_model():
multi = MultiProcessDLM(cp6, E2, models, order,
prior_model_prob,
- mean_prior=mean_prior,
- var_prior=var_prior,
+ m0=m0, C0=C0, n0=n0, s0=s0,
approx_steps=1)
return multi
@@ -56,15 +56,14 @@ def get_multi_model():
y = datasets.table_22()
x = [[1]]
- mean_prior = (0, 1)
- var_prior = (1, 0.01)
+ m0, C0 = (0, 1)
+ n0, s0 = (1, 0.01)
discounts = np.arange(0.7, 1.01, 0.1)
models = {}
for delta in discounts:
- models['%.2f' % delta] = DLM(y, x, mean_prior=mean_prior,
- var_prior=var_prior,
- discount=delta)
+ models['%.2f' % delta] = DLM(y, x, m0=m0, C0=C0, n0=n0, s0=s0,
+ state_discount=delta)
mix = DLMMixture(models)
multi = get_multi_model()
View
117 statlib/arma.py
@@ -2,13 +2,14 @@
from numpy.linalg import inv
import numpy as np
import numpy.linalg as LA
+import scipy.stats as stats
import statlib.plotting as plotting
+import statlib.dlm as dlm
-from pandas.util.testing import set_trace as st
+from pandas.util.testing import set_trace, debug
-from scikits.statsmodels.tsa.arima import ARMA as sm_arma
-from scikits.statsmodels.tsa.stattools import acf
from scikits.statsmodels.tools.decorators import cache_readonly
+import scikits.statsmodels.tsa.api as tsa_api
import scikits.statsmodels.api as sm
@@ -19,21 +20,40 @@ def __init__(self, data, p=1):
self.p = p
self.arx = data - data.mean()
- self.exog, self.endog = _prep_arvars(self.arx, p)
+ self.endog, self.exog = _prep_arvars(self.arx, p)
@cache_readonly
def ref_results(self):
- return sm.OLS(self.exog, self.endog).fit()
+ return tsa_api.AR(self.arx).fit(self.p, trend='nc')
@property
def dlm_rep(self):
"""
DLM form of AR(p) state transition matrix
"""
- beta = LA.lstsq(self.endog, self.exog)[0]
+ beta = LA.lstsq(self.exog, self.endog)[0]
lower = np.c_[np.eye(self.p-1), np.zeros((self.p-1, 1))]
return np.vstack((beta, lower))
+ def fit_dlm(self, m0, M0, n0, s0, discount=1.):
+ """
+ Fit dynamic linear model to the data, with the AR parameter as the state
+ vector
+
+ Parameters
+ ----------
+ m0 : ndarray (p)
+ M0 : ndarray (p x p)
+ n0 : int / float
+ s0 : float
+
+ Returns
+ -------
+ model : DLM
+ """
+ return dlm.DLM(self.endog, self.exog, mean_prior=(m0, M0),
+ var_prior=(n0, s0), discount=discount)
+
def decomp(self):
"""
Compute decomposition of input time series into components
@@ -53,10 +73,11 @@ def decomp(self):
decomp : n x k
"""
p = self.p
- X = self.endog
+ X = self.exog
vals, vecs = LA.eig(self.dlm_rep)
+ # (B_t'F) B_t^-1, where F = [1, 0, ..., 0]
H = np.dot(np.diag(vecs[0]), inv(vecs))
mods = np.abs(vals)
@@ -155,6 +176,78 @@ def ar_simulate(phi, s, n, dist='normal'):
return out
+def ar_model_like(data, maxp, m0, M0, n0, s0):
+ """
+ Compute approximate marginal likelihood for univariate AR(p) model for
+ p=0,...maxp.
+
+ Parameters
+ ----------
+ data
+ maxp : int
+ Maximum model order
+ m0 : ndarray (maxp x maxp)
+ Prior mean for state
+ M0 : ndarray (maxp x maxp)
+ Prior state variance matrix
+ n0 : int / float
+ Prior degrees of freedom for obs var
+ s0 : float
+ Prior estimate of observation variance
+
+ Notes
+ -----
+ p(\theta) ~ N(m0, M0)
+ \nu ~ IG(n0 / 2, n0 * s0 / 2)
+ """
+ from scipy.special import gammaln as gamln
+
+ arx = data - data.mean()
+ nobs = len(arx)
+
+ result = np.zeros(maxp + 1)
+ for p in xrange(maxp + 1):
+ llik = 0.
+ model = ARModel(data, p=p)
+
+ mt = m0[:p]; Mt = M0[:p, :p]
+ st = s0; nt = n0
+
+ # dlm_model = model.fit_dlm(m0[:p], M0[:p, :p], n0, s0)
+ # F = dlm_model.F
+ # y = dlm_model.y
+
+ # for t in xrange(p, nobs - p):
+ # nt = dlm_model.df[t]
+ # q = dlm_model.forc_var[t-1]
+ # st = dlm_model.var_est[t - 1]
+ # e = y[t - 1] - np.dot(F[t - 1], dlm_model.mu_mode[t - 1])
+
+ # # Student-t log pdf
+ # sd = np.sqrt(q * st)
+ # llik += stats.t(nt).logpdf(e / sd) - np.log(sd)
+
+ for t in xrange(p, nobs):
+ # DLM update equations
+ x = arx[t-p:t][::-1]
+ A = np.dot(Mt, x); q = np.dot(x, A) + 1; A /= q
+ e = arx[t] - np.dot(x, mt)
+
+ if t >= 2 * maxp:
+ sd = np.sqrt(q * st)
+ llik += stats.t(nt).logpdf(e / sd) - np.log(sd)
+ mt = mt + A * e
+ st = (nt * st + np.dot(e, e) / q) / (nt + 1)
+ nt = nt + 1
+ Mt = Mt - np.outer(A, A) * q
+
+ result[p] = llik
+
+ popt = result.argmax()
+ maxlik = result[popt]
+
+ return result
+
if __name__ == '__main__':
from statlib.linmod import BayesLM
import statlib.datasets as ds
@@ -168,7 +261,13 @@ def ar_simulate(phi, s, n, dist='normal'):
p = 8
y, X = _prep_arvars(eeg, p)
- beta_prior = (np.zeros(p), np.eye(p) * 100)
- var_prior = (2, 1000)
+ m0, M0 = np.zeros(p), np.eye(p) * 100
+ n0, s0 = 2, 1000
+ beta_prior = (m0, M0)
+ var_prior = (n0, s0)
model2 = BayesLM(y, X, beta_prior, var_prior)
+
+ maxp = 20
+ result = ar_model_like(eeg, maxp, np.zeros(maxp),
+ np.eye(maxp) * 1, 1, 1)
View
11 statlib/datasets.py
@@ -80,11 +80,20 @@ def foo():
return pn.DataFrame(datad)
-def eeg_data():
+def eeg400_data():
+ """
+ Cz series, 400 data points
+ """
data = open(os.path.join(data_path, 'eeg.dat')).read().strip()
data = re.sub('[\s\n]+', ' ', data).split()
return np.array(data, dtype=float)
+def eeg_data():
+ """
+ Cz series, 3600 data points
+ """
+ return np.loadtxt('statlib/data/eeg2.dat')
+
def parse_table_22():
path = os.path.join(data_path, 'Table2.2.data.txt')
sep = '\s+'
View
293 statlib/dlm.py
@@ -26,6 +26,8 @@
from statlib.tools import chain_dot, nan_array
import statlib.distributions as distm
+from pandas.util.testing import set_trace as st
+
m_ = np.array
def nct_pdf(x, df, nc):
@@ -34,25 +36,50 @@ def nct_pdf(x, df, nc):
return dt(x, df, nc)
class DLM(object):
- """
- Bayesian Gaussian Dynamic Linear Model (DLM) with discount factor
+ r"""
+ Bayesian Gaussian Dynamic Linear Model (DLM) with discount factors
Parameters
----------
y : ndarray n x 1
Response variable
- F : ndarray n x k
+ F : ndarray (n x k or 1d for constant Ft = F)
Regressor matrix
G : ndarray k x k
State transition matrix
- mean_prior : tuple (mean, variance)
- mean: length k
- variance: k x k
- Normal prior for mean response
- var_prior : tuple (a, b)
- Inverse-gamma prior for observation variance
- discount : float
+ m0 : ndarray k
+ Prior mean for state vector
+ C0 : ndarray k x k
+ Prior covariance for state vector
+ n0 : int
+ Prior degrees of freedom for observation variance
+ s0 : float
+ Prior point estimate of observation variance
+ state_discount : float
+ Discount factor for determining state evolution variance
Wt = Ct * (1 - d) / d
+ var_discount : float
+ Discount factor for observational variance
+
+ Notes
+ -----
+ Normal dynamic linear model (DLM) in state space form
+
+ .. math::
+
+ y = F_t' \theta_t + \nu_t \\
+ \theta_t = G_t \theta_{t-1} + \omega_t \\
+ \nu_t \sim N(0, V_t) \\
+ \omega_t \sim N(0, W_t) \\
+ W_t = C_t (1 - \delta) / \delta \\
+ IG(df0 / 2, df0 * v0 / 2)
+
+ Priors
+
+ .. math::
+
+ {\sf Normal}(m_0, C_0)\\
+ {\sf InverseGamma}(n_0 / 2, n_0 s_0 / 2)
"""
F = None
mu_mode = None # mt
@@ -62,47 +89,69 @@ class DLM(object):
var_est = None # St
forc_var = None # Qt
- def __init__(self, y, F, G=None, mean_prior=None, var_prior=None,
- discount=0.9):
- self.dates = y.index
- self.y = m_(y)
- self.nobs = len(y)
+ def __init__(self, y, F, G=None, m0=None, C0=None, n0=None, s0=None,
+ mean_prior=None, var_prior=None, discount=None,
+ state_discount=0.9, var_discount=1.):
+ try:
+ self.dates = y.index
+ except AttributeError:
+ self.dates = None
- if self.y.ndim == 1:
+ y = m_(y)
+ if y.ndim == 1:
pass
else:
raise Exception
+ self.y = y
+ self.nobs = len(y)
+
+ constant = False
+
F = m_(F)
if F.ndim == 1:
- self.ndim = len(F)
+ if len(F) == len(y):
+ self.ndim = 1
+ else:
+ constant = True
+ self.ndim = len(F)
else:
self.ndim = F.shape[1]
f_shape = self.nobs, self.ndim
+ # HACK ?
# constant DLM handling
if F.ndim == 1:
- self.F = np.ones(f_shape) * F
+ if constant:
+ F = np.ones(f_shape) * F
+ else:
+ F = F.reshape(f_shape)
else:
if len(F) == 1 and self.nobs > 1:
F = np.ones(f_shape) * F[0]
- self.F = F
+
+ self.F = F
if G is None:
G = np.eye(self.ndim)
self.G = G
- self.mean_prior = mean_prior
- self.var_prior = var_prior
- self.disc = discount
+ # HACK for transition
+ if discount is not None:
+ state_discount = discount
- df0, df0v0 = var_prior
- self.df = df0 + np.arange(self.nobs + 1) # precompute
+ self.state_discount = state_discount
+ self.var_discount = var_discount
- self.m0, self.C0 = mean_prior
- self.df0 = df0
- self.v0 = df0v0 / df0
+ if mean_prior is not None:
+ m0, C0 = mean_prior
+
+ if var_prior is not None:
+ n0, s0 = var_prior
+
+ self.m0, self.C0 = m0, C0
+ self.df0, self.s0 = n0, s0
self._compute_parameters()
@@ -124,11 +173,15 @@ def _compute_parameters(self):
(self.mu_mode,
self.mu_forc_mode,
self.mu_scale,
+ self.df,
self.var_est,
self.forc_var,
- self.mu_forc_scale) = _filter_python(self.y, self.F, self.G, self.disc,
- self.df0, self.v0, self.m0,
- self.C0)
+ self.mu_forc_scale) = _filter_python(self.y, self.F,
+ self.G,
+ self.state_discount,
+ self.var_discount,
+ self.df0, self.s0,
+ self.m0, self.C0)
def backward_sample(self, steps=1):
"""
@@ -224,9 +277,10 @@ def plot_mu(self, alpha=0.10, prior=True, index=None, ax=None):
@property
def forc_dist(self):
+ delta = self.state_discount
return stats.t(self.df[:-1],
loc=self.mu_mode[:-1],
- scale=self.mu_scale[:-1] / self.disc)
+ scale=self.mu_scale[:-1] / delta)
@property
def mupost_dist(self):
@@ -265,14 +319,14 @@ def mu_ci(self, alpha=0.10, prior=False):
diags = self.mu_scale[:, _x, _y]
# Only care about marginal scale
- delta = self.disc
+ delta = self.state_discount
if isinstance(delta, np.ndarray):
delta = np.diag(delta)
if prior:
df = self.df[:-1]
mode = self.mu_mode[:-1]
- scale = np.sqrt(diags[:-1] / self.disc)
+ scale = np.sqrt(diags[:-1] / delta)
else:
df = self.df[1:]
mode = self.mu_mode[1:]
@@ -289,19 +343,16 @@ def forc_ci(self, alpha=0.10):
return distm.make_t_ci(self.df[:-1], self.forecast,
np.sqrt(self.forc_var), alpha=alpha)
- def fit(self, discount=0.9):
- pass
+class VarDiscountDLM(DLM):
+ pass
class DLM2(DLM):
def _compute_parameters(self):
- df0, df0v0 = self.var_prior
- v0 = float(df0v0) / df0
-
- m0, C0 = self.mean_prior
-
(mode, a, Q, C, S) = _filter_cython(self.y, self.F, self.G,
- self.disc, df0, v0, m0, C0)
+ self.state_discount,
+ self.df0, self.s0,
+ self.m0, self.C0)
self.forc_var = Q
self.mu_forc_mode = a
@@ -457,7 +508,80 @@ def _mvfilter_python(Y, F, G, V, delta, beta, df0, v0, m0, C0):
return mode, a, C, S, Q
-def _filter_python(Y, F, G, delta, df0, v0, m0, C0):
+from line_profiler import LineProfiler
+prof = LineProfiler()
+
+class DLMFilter(object):
+ # should go into Cython
+ pass
+
+def _filter_python_multi(Y, F, G, delta, df0, v0, m0, C0):
+ """
+ Univariate DLM update equations with unknown observation variance
+
+ Done with numpy.matrix objects so works with multivariate case. About 3-4x
+ slower than specializing to the 1d case and using vectors / regular ndarrays
+ """
+ nobs = len(Y)
+ ndim = len(G)
+
+ mode = nan_array(nobs + 1, ndim)
+ a = nan_array(nobs + 1, ndim)
+ C = nan_array(nobs + 1, ndim, ndim)
+ S = nan_array(nobs + 1)
+ Q = nan_array(nobs)
+ R = nan_array(nobs + 1, ndim, ndim)
+
+ df = df0 + np.arange(nobs + 1) # precompute
+ S[0] = St = float(v0)
+
+ mode[0] = m0
+ C[0] = Ct = np.asmatrix(C0)
+
+ mt = np.asmatrix(m0).T
+ Y = np.asmatrix(Y).T
+ F = np.asmatrix(F)
+ G = np.asmatrix(G)
+
+ # allocate result arrays
+ for i in xrange(nobs):
+ # column vector, for W&H notational consistency
+ Yt = Y[i].T
+ Ft = F[i].T
+
+ # advance index: y_1 through y_nobs, 0 is prior
+ t = i + 1
+
+ # derive innovation variance from discount factor
+ at = mt
+ Rt = Ct
+ if t > 1:
+ # only discount after first time step!
+ if G is not None:
+ at = G * mt
+ Rt = G * Ct * G.T / delta
+ else:
+ Rt = Ct / delta
+
+ Q[t-1] = Qt = Ft.T * Rt * Ft + St
+ At = Rt * Ft / Qt
+
+ # forecast theta as time t
+ ft = Ft.T * at
+ e = Yt - ft
+
+ # update mean parameters
+ mt = at + At * e
+ S[t] = St = St + (St / df[t]) * (e * e.T / Qt - 1)
+ C[t] = Ct = (St / S[t-1]) * (Rt - At * At.T * Qt)
+
+ mode[t] = mt.T
+ a[t] = at.T
+ R[t] = Rt
+
+ return mode, a, C, S, Q, R
+
+def _filter_python_old(Y, F, G, delta, df0, v0, m0, C0):
"""
Univariate DLM update equations with unknown observation variance
"""
@@ -501,18 +625,87 @@ def _filter_python(Y, F, G, delta, df0, v0, m0, C0):
# forecast theta as time t
ft = np.dot(Ft.T, at)
- err = obs - ft
+ e = obs - ft
# update mean parameters
- mode[t] = at + np.dot(At, err)
- S[t] = S[t-1] + (S[t-1] / df[t]) * ((err ** 2) / Qt - 1)
+ mode[t] = at + np.dot(At, e)
+ S[t] = S[t-1] + (S[t-1] / df[t]) * ((e ** 2) / Qt - 1)
C[t] = (S[t] / S[t-1]) * (Rt - np.dot(At, At.T) * Qt)
a[t] = at
Q[t-1] = Qt
R[t] = Rt
- return mode, a, C, S, Q, R
+ return mode, a, C, df, S, Q, R
+
+def _filter_python(Y, F, G, delta, beta, df0, v0, m0, C0):
+ """
+ Univariate DLM update equations with unknown observation variance
+
+ delta : state discount
+ beta : variance discount
+ """
+ nobs = len(Y)
+ ndim = len(G)
+
+ mode = nan_array(nobs + 1, ndim)
+ a = nan_array(nobs + 1, ndim)
+ C = nan_array(nobs + 1, ndim, ndim)
+ S = nan_array(nobs + 1)
+ Q = nan_array(nobs)
+ R = nan_array(nobs + 1, ndim, ndim)
+ df = nan_array(nobs + 1)
+
+ mode[0] = mt = m0
+ C[0] = Ct = C0
+ df[0] = nt = float(df0)
+ S[0] = St = v0
+
+ dt = df0 * v0
+
+ # allocate result arrays
+ for i, obs in enumerate(Y):
+ # column vector, for W&H notational consistency
+ Ft = F[i:i+1].T
+
+ # advance index: y_1 through y_nobs, 0 is prior
+ t = i + 1
+
+ # derive innovation variance from discount factor
+ at = mt
+ Rt = Ct
+ if t > 1:
+ # only discount after first time step?
+ if G is not None:
+ at = np.dot(G, mt)
+ Rt = chain_dot(G, Ct, G.T) / delta
+ else:
+ Rt = Ct / delta
+
+ Qt = chain_dot(Ft.T, Rt, Ft) + St
+ At = np.dot(Rt, Ft) / Qt
+
+ # forecast theta as time t
+ ft = np.dot(Ft.T, at)
+ e = obs - ft
+
+ # update mean parameters
+ mode[t] = mt = at + np.dot(At, e)
+ dt = beta * dt + St * e * e / Qt
+ nt = beta * nt + 1
+ St = dt / nt
+
+ S[t] = St
+ Ct = (S[t] / S[t-1]) * (Rt - np.dot(At, At.T) * Qt)
+ Ct = (Ct + Ct.T) / 2 # symmetrize
+
+ df[t] = nt
+ C[t] = Ct
+ a[t] = at
+ Q[t-1] = Qt
+ R[t] = Rt
+
+ return mode, a, C, df, S, Q, R
def rmatnorm(M, U, V):
"""
@@ -577,8 +770,12 @@ def _filter_cython(Y, F, G, delta, df0, v0, m0, C0):
m0 = np.zeros((p, k))
m0[1:k+1] = np.eye(k)
- mean_prior = m0, 1000 * np.eye(p)
- var_prior = 2, D0
+ C0 = 1000 * np.eye(p)
+
+ n0 = 2
+
+ mean_prior = m0, C0
+ var_prior = n0, D0
model = MVDLM(Y, F,
state_discount=delta, cov_discount=beta,
View
17 statlib/multi.py
@@ -143,8 +143,8 @@ class MultiProcessDLM(object):
cf. W&H Section 12.4
"""
def __init__(self, y, F, models, order, prior_model_prob,
- mean_prior=None, var_prior=None, approx_steps=1.):
-
+ m0=None, C0=None, n0=None, s0=None, approx_steps=1.):
+ self.approx_steps = approx_steps
self.dates = y.index
self.y = np.array(y)
self.nobs = len(y)
@@ -171,9 +171,6 @@ def __init__(self, y, F, models, order, prior_model_prob,
self.approx_steps = 1
# self.approx_steps = int(approx_steps)
- self.mean_prior = mean_prior
- self.var_prior = var_prior
-
# set up result storage for all the models
self.marginal_prob = nan_array(self.nobs + 1, self.nmodels)
self.post_prob = nan_array(self.nobs + 1,
@@ -190,15 +187,15 @@ def __init__(self, y, F, models, order, prior_model_prob,
# set in initial values
self.marginal_prob[0] = self.prior_model_prob
- self.mu_mode[0], self.mu_scale[0] = mean_prior
+ self.mu_mode[0] = m0
+ self.mu_scale[0] = C0
# observation variance stuff
- n, d = var_prior
- self.df = n + np.arange(self.nobs + 1) # precompute
+ self.df = n0 + np.arange(self.nobs + 1) # precompute
self.var_est = nan_array(self.nobs + 1, self.nmodels)
self.var_scale = nan_array(self.nobs + 1, self.nmodels)
- self.var_est[0] = d / n
- self.var_scale[0] = d
+ self.var_est[0] = s0
+ self.var_scale[0] = s0 * n0
# forecasts are computed via mixture for now
self.forc_var = nan_array(self.nobs, self.nmodels, self.nmodels)

0 comments on commit e68ac66

Please sign in to comment.