Merge pull request #968 from springcoil/lightspeed

ENH: Timeseries GARCH and Lightspeed model

Author: twiecki

pymc3/examples/garch_example.py

@@ -1,9 +1,9 @@
 from pymc3 import Normal, sample, Model, Bound
-import theano.tensor as tt
+import theano.tensor as T
 import numpy as np
-import math
+
 """
-Example from STAN
+Example from STAN - slightly altered
 // GARCH(1,1)
 
 data {
@@ -29,27 +29,28 @@
   r ~ normal(mu,sigma);
 }
 """
-T = 8
-r = np.array([28,  8, -3,  7, -1,  1, 18, 12])
-sigma1 = np.array([15, 10, 16, 11,  9, 11, 10, 18])
-alpha0 = np.array([10, 10, 16, 8,  9, 11, 12, 18])
+J = 8
+r = np.array([28, 8, -3, 7, -1, 1, 18, 12])
+sigma1 = np.array([15, 10, 16, 11, 9, 11, 10, 18])
+alpha0 = np.array([10, 10, 16, 8, 9, 11, 12, 18])
 
 with Model() as garch:
-
-    alpha1 = Normal('alpha1', 0, 1, shape=T)
-    BoundedNormal = Bound(Normal, upper=(1-alpha1))
+    alpha1 = Normal('alpha1', 0, 1, shape=J)
+    BoundedNormal = Bound(Normal, upper=(1 - alpha1))
     beta1 = BoundedNormal('beta1', 0, sd=1e6)
     mu = Normal('mu', 0, sd=1e6)
 
-    theta = tt.sqrt(alpha0 + alpha1*tt.pow(r - mu, 2) + beta1 * tt.pow(sigma1, 2))
+    theta = T.sqrt(alpha0 + alpha1 * T.pow(r - mu, 2) + beta1 * T.pow(sigma1, 2))
 
     obs = Normal('obs', mu, sd=theta, observed=r)
 
+
 def run(n=1000):
     if n == "short":
         n = 50
     with garch:
         tr = sample(n)
 
+
 if __name__ == '__main__':
-    run()
+    run()

pymc3/examples/lightspeed_example.py

@@ -0,0 +1,44 @@
+import numpy as np
+import pymc3 as pm
+
+light_speed = np.array([28, 26, 33, 24, 34, -44, 27, 16, 40, -2, 29, 22, 24, 21, 25,
+                        30, 23, 29, 31, 19, 24, 20, 36, 32, 36, 28, 25, 21, 28, 29,
+                        37, 25, 28, 26, 30, 32, 36, 26, 30, 22, 36, 23, 27, 27, 28,
+                        27, 31, 27, 26, 33, 26, 32, 32, 24, 39, 28, 24, 25, 32, 25,
+                        29, 27, 28, 29, 16, 23])
+
+model_1 = pm.Model()
+
+with model_1:
+    # priors as specified in stan model
+    # mu = pm.Uniform('mu', lower = -tt.inf, upper= np.inf)
+    # sigma = pm.Uniform('sigma', lower = 0, upper= np.inf)
+
+
+    # using vague priors works
+    mu = pm.Uniform('mu', lower=light_speed.std() / 1000.0, upper=light_speed.std() * 1000.0)
+    sigma = pm.Uniform('sigma', lower=light_speed.std() / 1000.0, upper=light_speed.std() * 1000.0)
+
+    # define likelihood
+    y_obs = pm.Normal('Y_obs', mu=mu, sd=sigma, observed=light_speed)
+
+    # inference fitting the model
+
+    # I have to use slice because the following command
+    # trace = pm.sample(5000)
+    # produce the error
+    # ValueError: Cannot construct a ufunc with more than 32 operands
+    # (requested number were: inputs = 51 and outputs = 1)valueerror
+
+
+def run(n=5000):
+    with model_1:
+        xstart = pm.find_MAP()
+        xstep = pm.Slice()
+        trace = pm.sample(5000, xstep, xstart, random_seed=123, progressbar=True)
+
+        pm.summary(trace)
+
+
+if __name__ == '__main__':
+    run()