PySSM: Bayesian Inference of Linear Gaussian State Space Models

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

Authors:	Christopher Strickland, Robert Burdett, Kerrie Mengersen, Robert Denham
Title:	PySSM: A Python Module for Bayesian Inference of Linear Gaussian State Space Models
Abstract:	PySSM is a Python package that has been developed for the analysis of time series using linear Gaussian state space models. PySSM is easy to use; models can be set up quickly and efficiently and a variety of different settings are available to the user. It also takes advantage of scientific libraries NumPy and SciPy and other high level features of the Python language. PySSM is also used as a platform for interfacing between optimized and parallelized Fortran routines. These Fortran routines heavily utilize basic linear algebra and linear algebra Package functions for maximum performance. PySSM contains classes for filtering, classical smoothing as well as simulation smoothing.
	Page views:: 5207. Submitted: 2012-03-23. Published: 2014-04-07.
Paper:	PySSM: A Python Module for Bayesian Inference of Linear Gaussian State Space Models     [size=1em]Download PDF (Downloads: 5736)
Supplements:	pyssm-1.1e.tar.gz: Python source package [size=1em]Download(Downloads: 277; 95KB) pyssm-1.1e.win32-py2.7.exe: Windows 32-bit binary package [size=1em]Download(Downloads: 254; 2MB) v57i06-examples.zip: Python example code from the paper [size=1em]Download(Downloads: 270; 5KB)

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏1 回帖

关键词：State Space Inference Bayesian Gaussian models scientific different advantage available developed

1. from numpy import log, ones, column_stack, hstack, mean
   
2. from numpy import random, zeros, sqrt
   
3. from pymcmc.mcmc import MCMC, CFsampler, IndMH
   
4. from pymcmc.regtools import LinearModel, CondScaleSampler
   
5. import pylab
   
6. from pyssm.ssm import Filter, SimSmoother
   
7. 
   
8. 
   
9. def simdata(nobs):
   
10.     ht = 1.0 ** 2
    
11.     zt = 1.0
    
12.     tt = 0.95
    
13.     qt = 0.3 ** 2
    
14.     gt = 1.0
    
15.     rt = 1.0
    
16.     mu = 5.0
    
17.     a1 = mu
    
18.     p1 = qt / (1. - tt ** 2)
    
19.     beta = mu * (1.0 - tt)
    
20.     wmat = ones((1, 1, nobs))
    
21. 
    
22.     filt = Filter(nobs, 1, 1, 1, False, wmat=wmat)
    
23.     filt.initialise_system(a1, p1, zt, ht, tt, gt, qt, rt, beta=beta)
    
24.     filt.simssm()
    
25. 
    
26.     yvec = filt.get_ymat().T.flatten()
    
27.     simstate = filt.get_state()
    
28.     return yvec, simstate
    
29. 
    
30. 
    
31. def simstate(store):
    
32.     system = store['simsm'].get_system()
    
33.     system.update_tt(store['rho'])
    
34. 
    
35.     p1 = system._qt() / (1.0 - store['rho'] ** 2)
    
36.     a1 = store['beta'] / (1.0 - store['rho'])
    
37. 
    
38.     system.update_a1(a1)
    
39.     system.update_p1(p1)
    
40. 
    
41.     system.update_beta(store['beta'])
    
42. 
    
43.     state = store['simsm'].sim_smoother()
    
44. 
    
45.     return state
    
46. 
    
47. 
    
48. def simht(store):
    
49.     system = store['simsm'].get_system()
    
50.     residual = store['simsm'].get_meas_residual()
    
51.     ht = store['scale_sampler'].sample(residual.T)
    
52.     system.update_ht(ht ** 2)
    
53. 
    
54.     return ht ** 2
    
55. 
    
56. 
    
57. def prior_rho(store):
    
58.     if store['rho'] < 1.0 and store['rho'] > 0.0:
    
59.         rho1 = 15.0
    
60.         rho2 = 1.5
    
61.         return (rho1 - 1.0) * log(store['rho']) + (rho2 - 1.0) * log(1.0 - store['rho'])
    
62.     else:
    
63.         return -1E256
    
64. 
    
65. 
    
66. def prior_sigma(store):
    
67.     nu = 10.0
    
68.     S = 0.01
    
69.     return -(nu + 1) * log(store['sigma']) - S / (2.0 * store['sigma'] ** 2)
    
70. 
    
71. 
    
72. def post_rho_sigma(store):
    
73.     if store['rho'] < 1.0 and store['rho'] > 0.0:
    
74.         system = store['simsm'].get_system()
    
75.         p1 = system._qt() / (1.0 - store['rho'] ** 2)
    
76.         a1 = store['beta'] / (1.0 - store['rho'])
    
77.         system.update_p1(p1)
    
78.         system.update_a1(a1)
    
79.         system.update_tt(store['rho'])
    
80.         system.update_beta(store['beta'])
    
81.         system.update_qt(store['sigma'] ** 2)
    
82.         lnpr = store['simsm'].log_probability_state()
    
83. 
    
84.         #lnpr =  store['simsm'].calc_log_likelihood()
    
85.         return lnpr + prior_rho(store) + prior_sigma(store)
    
86.     else:
    
87.         return -1E256
    
88. 
    
89. 
    
90. def cand_rho_sigma(store):
    
91.     nobs = store['state'].shape[1]
    
92.     store['bayes_reg'].update_yvec(store['state'][0, 1:nobs])
    
93.     xmat = column_stack([ones(nobs - 1), store['state'][0, 0:nobs - 1]])
    
94.     store['bayes_reg'].update_xmat(xmat)
    
95.     sig, beta = store['bayes_reg'].sample()
    
96.     return sig, beta[0], beta[1]
    
97. 
    
98. 
    
99. def cand_prob(store):
    
100.     beta = hstack([store['beta'], store['rho']])
     
101.     return store['bayes_reg'].log_posterior_probability(store['sigma'], beta)
     
102. 
     
103. 
     
104. def transform_beta(store):
     
105.     return store['beta'] / (1.0 - store['rho'])
     
106. 
     
107. 
     
108. def main():
     
109.     """
     
110.     The main routine.
     
111.     """
     
112. 
     
113.     random.seed(12345)
     
114.     nobs = 1000
     
115.     nstate = 1
     
116.     rstate = 1
     
117.     yvec, simulated_state = simdata(nobs)
     
118.     data = {'yvec': yvec}
     
119. 
     
120.     ht = 1.0
     
121.     tt = 0.9
     
122.     zt = 1.0
     
123.     qt = 0.4
     
124.     gt = 1.0
     
125.     rt = 1.0
     
126.     p1 = qt / (1. - tt ** 2)
     
127.     mu = mean(yvec)
     
128.     a1 = mu
     
129.     beta = mu * (1 - tt)
     
130.     wmat = 1.0
     
131. 
     
132.     data['simsm'] = SimSmoother(yvec, nstate, rstate,
     
133.                                 False, properties={'gt': 'eye'},
     
134.                                 wmat=wmat)
     
135.     data['simsm'].initialise_system(a1, p1, zt, ht, tt, gt, qt, rt, beta=beta)
     
136.     data['scale_sampler'] = CondScaleSampler(prior=['inverted_gamma', 10.0, 0.01])
     
137.     data['bayes_reg'] = LinearModel(yvec[1:nobs], column_stack([ones(nobs - 1), yvec[0:nobs - 1]]))
     
138.     samplestate = CFsampler(simstate, zeros((nstate, nobs)), 'state')
     
139.     sampleht = CFsampler(simht, ht, 'ht')
     
140.     samplesigbetarho = IndMH(cand_rho_sigma, post_rho_sigma,
     
141.                              cand_prob, [sqrt(qt), 0.1, tt],
     
142.                              ['sigma', 'beta', 'rho'])
     
143. 
     
144.     blocks1 = [samplestate, sampleht, samplesigbetarho]
     
145.     mcmc = MCMC(5000, 2000, data, blocks1,
     
146.                 transform={'beta': transform_beta}, runtime_output = True)
     
147.     mcmc.sampler()
     
148.     mcmc.output(parameters=['ht', 'sigma', 'rho', 'beta'],
     
149.                filename='ar1out.txt')
     
150. 
     
151.     means, vars = mcmc.get_mean_var('state')
     
152.     pylab.plot(range(nobs), means[0], color='k')
     
153.     pylab.plot(range(nobs), simulated_state[0], color='k',)
     
154.     pylab.title('Simulated vs estimated state')
     
155.     pylab.savefig('AR1.pdf')
     
156.     pylab.close("all")
     
157.    
     
158. if __name__ == '__main__':
     
159.     main()

复制代码

1. import numpy as np
   
2. import os
   
3. import pymcmc.mcmc as mcmc
   
4. import pymcmc.regtools as regtools
   
5. import pylab
   
6. import pyssm.ssm as ssm
   
7. 
   
8. """ Get the path for the data. If this was
   
9. installed using setup.py it will be in the
   
10. data directory of the module"""
    
11. 
    
12. datadir = os.path.join(os.path.dirname(ssm.__file__), 'data')
    
13. 
    
14. 
    
15. def update_tt(store):
    
16.     system = store['simsm'].get_system()
    
17.     tt = system.tt()
    
18.     tt[2, 2] = store['rho'] * np.cos(store['lambda'])
    
19.     tt[2, 3] = store['rho'] * np.sin(store['lambda'])
    
20.     tt[3, 2] = -tt[2, 3]
    
21.     tt[3, 3] = tt[2, 2]
    
22.     system.update_tt(tt)
    
23. 
    
24. 
    
25. def update_p1(store):
    
26.     system = store['simsm'].get_system()
    
27.     p1 = system.p1()
    
28.     p1[2, 2] = store['sigma_cycle'] ** 2 / \
    
29.                (1. - store['rho'] ** 2)
    
30.     p1[3, 3] = p1[2, 2]
    
31.     system.update_p1(p1)
    
32. 
    
33. 
    
34. def update_qt(store):
    
35.     system = store['simsm'].get_system()
    
36.     qt = np.zeros((3, 3))
    
37.     qt[0, 0] = store['sigma_level'] ** 2
    
38.     qt[1, 1] = store['sigma_cycle'] ** 2
    
39.     qt[2, 2] = qt[1, 1]
    
40.     system.update_qt(qt)
    
41. 
    
42. 
    
43. def simstate(store):
    
44.     update_tt(store)
    
45.     update_qt(store)
    
46.     update_p1(store)
    
47.     return store['simsm'].sim_smoother()
    
48. 
    
49. 
    
50. def simht(store):
    
51.     residual = store['simsm'].get_meas_residual()
    
52.     ht = np.linalg.inv(
    
53.         store['scale_sampler'].sample(residual.T)
    
54.         )
    
55.     system = store['simsm'].get_system()
    
56.     system.update_ht(ht)
    
57.    
    
58.     return ht
    
59. 
    
60. 
    
61. def simsig_level(store):
    
62.     update_tt(store)
    
63.     residual = store['simsm'].get_state_residual(
    
64.         state_index=[0])
    
65.     sigma_level = store['scale_sampler2'].sample(
    
66.         residual.T)
    
67.     return sigma_level
    
68. 
    
69. 
    
70. def posterior_sig_cycle(store):
    
71.     if store['sigma_cycle'] > 0:
    
72.         update_tt(store)
    
73.         update_qt(store)
    
74.         update_p1(store)
    
75.         probstate = store['simsm'].log_probability_state()
    
76.         #probstate = store['simsm'].log_likelihood()
    
77.         return probstate + prior_sig(store['sigma_cycle'])
    
78.     else:
    
79.         return -1E256
    
80. 
    
81. 
    
82. def prior_sig(sig):
    
83.     nu = -1.
    
84.     S = 0.0
    
85.     return -(nu + 1) * np.log(sig) - S / (2.0 * sig ** 2)
    
86. 
    
87. 
    
88. def posterior_lambda(store):
    
89.     update_tt(store)
    
90.     lnpr = store['simsm'].log_probability_state()
    
91.     #lnpr = store['simsm'].log_likelihood()
    
92.     return lnpr + prior_lambda(store)
    
93. 
    
94. 
    
95. def prior_lambda(store):
    
96.     """flat prior restricting the period of the cycle to be between 10
    
97.     months and 14 months"""
    
98.     if store['lambda'] > np.pi / 20. and store['lambda'] < 2 * np.pi / 4.:
    
99.         return 0.0
    
100.     else:
     
101.         return -1E256
     
102. 
     
103. 
     
104. def posterior_rho(store):
     
105.     if store['rho'] > 0 and store['rho'] < 1.0:
     
106.         update_tt(store)
     
107.         update_p1(store)
     
108.         lnpr = store['simsm'].log_probability_state()
     
109.         #lnpr = store['simsm'].log_likelihood()
     
110.         return lnpr + prior_rho(store)
     
111.     else:
     
112.         return -1E256
     
113. 
     
114. 
     
115. def prior_rho(store):
     
116.     rho1 = 15.0
     
117.     rho2 = 1.5
     
118.     return (rho1 - 1.0) * np.log(store['rho']) + (rho2 - 1.) * \
     
119.            np.log(1. - store['rho'])
     
120. 
     
121. 
     
122. def main():
     
123.     np.random.seed(12345)
     
124.     filename = os.path.join(datadir, 'farmb.txt')
     
125.     ymat = np.loadtxt(filename).T / 1000.
     
126.     nseries, nobs = ymat.shape
     
127.     data = {'ymat': ymat}
     
128. 
     
129.     nstate = 4
     
130.     rstate = 3
     
131.     ht = np.eye(nseries)
     
132.     zt = np.column_stack([np.ones(nseries), np.zeros(nseries),
     
133.                           np.ones(nseries), np.zeros(nseries)])
     
134.     rt = np.ones(nseries)
     
135.     tt = np.zeros((nstate, nstate))
     
136.     rho = 0.9
     
137.     lamb = 2. * np.pi / 20.
     
138.     tt[0, 0] = 1.0
     
139.     tt[0, 1] = 1.0
     
140.     tt[1, 1] = 1.0
     
141.     tt[2, 2] = rho * np.cos(lamb)
     
142.     tt[2, 3] = rho * np.sin(lamb)
     
143.     tt[3, 2] = -tt[2, 3]
     
144.     tt[3, 3] = tt[2, 2]
     
145. 
     
146.     sig_cycle = 0.3
     
147.     sig_level = 0.3
     
148.     qt = np.diag(np.array([sig_level, sig_cycle, sig_cycle]) ** 2)
     
149.     gt = np.zeros((nstate, rstate))
     
150.     gt[0, 0] = 1.0
     
151.     gt[2:, 1:] = np.eye(2)
     
152.     a1 = np.array([7.5, 0.000, 0.0, 0.0])
     
153.     p1 = np.zeros((nstate, nstate))
     
154.     p1[0, 0] = 10.
     
155.     p1[1, 1] = 10.
     
156.     p1[2, 2] = qt[2, 2] / (1 - rho ** 2)
     
157.     p1[3, 3] = p1[2, 2]
     
158. 
     
159.     timevar = False
     
160.     data['simsm'] = ssm.SimSmoother(ymat, nstate, rstate, timevar,
     
161.                                 properties={'rt': 'eye'})
     
162.     data['simsm'].initialise_system(a1, p1, zt, ht, tt, gt, qt, rt)
     
163. 
     
164.     prior_wishart = ['wishart', 10 * np.ones(nseries), 0.1 * np.eye(nseries)]
     
165.     data['scale_sampler'] = regtools.CondScaleSampler(prior=prior_wishart)
     
166.     data['scale_sampler2'] = regtools.CondScaleSampler(
     
167.         prior=['inverted_gamma', 10, 0.1])
     
168. 
     
169.     samplestate = mcmc.CFsampler(simstate,
     
170.                                  np.zeros((nstate, nobs)),
     
171.                                  'state', store='none')
     
172.     sampleht = mcmc.CFsampler(simht, ht, 'ht')
     
173.     samplesig_level = mcmc.CFsampler(simsig_level,
     
174.                                      sig_level, 'sigma_level')
     
175.     samplesig_cycle = mcmc.RWMH(posterior_sig_cycle, 0.05,
     
176.                                 sig_cycle ** 2, 'sigma_cycle',
     
177.                                 adaptive='GFS')
     
178.     samplerho = mcmc.RWMH(posterior_rho, 0.09, rho, 'rho', adaptive='GFS')
     
179.     samplelambda = mcmc.RWMH(posterior_lambda, 1.03, lamb,
     
180.                              'lambda', adaptive='GFS')
     
181. 
     
182.     blocks = [samplestate, sampleht, samplesig_cycle,
     
183.               samplesig_level, samplerho, samplelambda]
     
184.     mcmcobj = mcmc.MCMC(8000, 3000, data, blocks, runtime_output=True)
     
185.     mcmcobj.sampler()
     
186.     mcmcobj.output(parameters=['rho', 'lambda', 'sigma_level', 'sigma_cycle'])
     
187. 
     
188.     means, vars = mcmcobj.get_mean_var('state')
     
189. 
     
190.     pylab.subplot(2, 1, 1)
     
191.     pylab.title("Trend Cycle Model")
     
192.     pylab.plot(ymat.T, color='k')
     
193.     pylab.plot(means[0], color='k')
     
194.     pylab.subplot(2, 1, 2)
     
195.     pylab.plot(means[2], color='k')
     
196.     pylab.savefig('trendcycle.pdf')
     
197.     pylab.close("all")
     
198.    
     
199. if __name__ == '__main__':
     
200.     main()

复制代码

1. import numpy as np
   
2. import os
   
3. import pylab
   
4. 
   
5. import pymcmc.mcmc as mcmc
   
6. import pymcmc.regtools as regtools
   
7. 
   
8. import pyssm.ssm as ssm
   
9. 
   
10. """ Get the path for the data. If this was installed using setup.py
    
11. it will be in the data directory of the module"""
    
12. datadir = os.path.join(os.path.dirname(ssm.__file__), 'data')
    
13. 
    
14. 
    
15. def update_gt(store):
    
16.     system = store['simsm'].get_system()
    
17.     gt = np.eye(2) * store['sigma']
    
18.     system.update_gt(gt)
    
19. 
    
20. 
    
21. def simstate(store):
    
22.     update_gt(store)
    
23.     #import pdb
    
24.     #pdb.set_trace()
    
25.     state = store['simsm'].sim_smoother()
    
26.     return state
    
27. 
    
28. 
    
29. def sim_ht(store):
    
30.     residual = store['simsm'].get_meas_residual()
    
31.     sqrtht = store['scale_sampler'].sample(residual.T)
    
32.     system = store['simsm'].get_system()
    
33.     ht = sqrtht ** 2
    
34.     system.update_ht(ht)
    
35.     return ht
    
36. 
    
37. 
    
38. def post_sigma(store):
    
39.     if store['sigma'] > 0:
    
40.         update_gt(store)
    
41.         lnpr = store['simsm'].log_likelihood()
    
42.         #lnpr = store['simsm'].log_probability_state(diffuse = True)
    
43. 
    
44.         return lnpr + prior_sigma(store['sigma'])
    
45.     else:
    
46.         return -1E256
    
47. 
    
48. 
    
49. def prior_sigma(sigma):
    
50.     nu = 10
    
51.     S = 0.1
    
52.     return -(nu + 1) * np.log(sigma) - 0.5 * S / sigma ** 2
    
53. 
    
54. 
    
55. def main():
    
56.     np.random.seed(1234)
    
57.     data = np.loadtxt(os.path.join(datadir, 'motorcycle.txt'))
    
58.     data = data[data[:, 1] > 0]
    
59.     yvec = data[:, 2]
    
60.     delta = data[:, 1]
    
61.     nobs = yvec.shape[0]
    
62.     nstate = 2
    
63.     rstate = 2
    
64.     data = {'yvec': yvec}
    
65. 
    
66.     sigeps = 0.9
    
67.     ht = sigeps ** 2
    
68.     rt = 1.
    
69.     zt = np.array([[1.0, 0.0]])
    
70.     tt = np.zeros((2, 2, nobs))
    
71.     tt[0, 0, :] = 1.
    
72.     tt[1, 1, :] = 1.
    
73.     tt[0, 1, :] = delta
    
74. 
    
75.     sigma = 0.3
    
76.     qt = np.zeros((2, 2, nobs))
    
77.     qt[0, 0, :] = delta ** 3 / 3.
    
78.     qt[0, 1, :] = qt[1, 0, :] = delta ** 2 / 2
    
79.     qt[1, 1, :] = delta
    
80. 
    
81.     gt = np.eye(nstate) * sigma
    
82.     a1 = np.zeros(2)
    
83.     p1 = np.eye(2)
    
84.     wmat = np.zeros((2, 2, nobs))
    
85. 
    
86.     data['simsm'] = ssm.SimSmoother(yvec, nstate, rstate,
    
87.                            timevar={'tt': True, 'qt': True},
    
88.                            properties={'rt': 'eye'}, wmat=wmat,
    
89.                            joint_sample=['diffuse', np.eye(2)])
    
90. 
    
91.     data['simsm'].initialise_system(a1, p1, zt, ht, tt, gt, qt, rt)
    
92.     data['scale_sampler'] = regtools.CondScaleSampler(
    
93.         prior=['inverted_gamma', 10, 0.1]
    
94.         )
    
95. 
    
96.     samplestate = mcmc.CFsampler(simstate, np.zeros((nstate, nobs)), 'state',
    
97.                                  store='none')
    
98.     sample_ht = mcmc.CFsampler(sim_ht, np.var(yvec), 'ht')
    
99.     sample_sigma = mcmc.RWMH(post_sigma, 1., 0.2, 'sigma', adaptive='GFS')
    
100. 
     
101.     blocks = [samplestate, sample_ht, sample_sigma]
     
102.     mcmcobj = mcmc.MCMC(8000, 3000, data, blocks, runtime_output=True)
     
103.     mcmcobj.sampler()
     
104.     mcmcobj.output(parameters=['ht', 'sigma'], filename='spline.out')
     
105. 
     
106.     means, vars = mcmcobj.get_mean_var('state')
     
107.     pylab.plot(np.cumsum(delta), yvec, '.', color='k')
     
108.     pylab.title("Smoothing spline for motorcycle acceleration data.")
     
109.     pylab.plot(np.cumsum(delta), means[0], color='k')
     
110.     pylab.savefig('spline.pdf')
     
111.     pylab.close("all")
     
112. 
     
113. if __name__ == '__main__':
     
114.     main()

复制代码

1. #!/usr/bin/env python
   
2. import os
   
3. 
   
4. try:
   
5.     os.system("clear")
   
6. except:
   
7.     os.system("cls")
   
8. 
   
9. 
   
10. intro = """
    
11. 
    
12. PySSM Examples
    
13. ==============
    
14. 
    
15. This script will run the three examples presented in the
    
16. accompanying paper. Scripts can also be run individually, using, for
    
17. example,
    
18. 
    
19. python example_ssm_ar1_reg.py
    
20. 
    
21. """
    
22. 
    
23. print intro
    
24. 
    
25. print """
    
26. 
    
27. Example 1: Autoregressive model with regressors
    
28. (see Section 4.1).
    
29. 
    
30. This script will produce a pdf called AR1.pdf (see Figure 1 in the paper),
    
31. and a summary of the model fit and parameters in a text file called ar1out.txt.
    
32. 
    
33. """
    
34. 
    
35. os.system("python example_ssm_ar1_reg.py")
    
36. 
    
37. print """
    
38. 
    
39. Example 2: Spline smoothing
    
40. (see Section 4.2)
    
41. 
    
42. This script will produce a pdf called spline.pdf (see Figure 2 in the paper),
    
43. and a summary of the model fit and parameters in a text file called spline.out
    
44. 
    
45. """
    
46. 
    
47. os.system("python example_ssm_vector_spline.py")
    
48. 
    
49. 
    
50. print """
    
51. 
    
52. Example 3: Trend plus cycle model
    
53. (see Section 4.3)
    
54. 
    
55. This script will produce a pdf called trendcycle.pdf  (see Figure 3 in the paper),
    
56. and a summary of the model fit and parameters will be printed to the screen.
    
57. 
    
58. """
    
59. 
    
60. os.system("python example_ssm_trendcycle2.py")

复制代码