title: "Bracket nesting level exceeded maximum" when working with large dataset in PyMC3
body: <p>I have a recursive model in PyMC3 and I'm training it on about 600 timesteps. I'm getting the error</p>
  <pre><code>Exception: ('Compilation failed (return status=1): /Users/tiwalayoaina/.theano/compiledir_macOS-10.14.6-x86_64-i386-64bit-i386-3.8.8-64/tmp1dgcmqm3/mod.cpp:30412:32: fatal error: bracket nesting level exceeded maximum of 256.         if (!PyErr_Occurred()) {.                                ^. /Users/actinidia/.theano/compiledir_macOS-10.14.6-x86_64-i386-64bit-i386-3.8.8-64/tmp1dgcmqm3/mod.cpp:30412:32: note: use -fbracket-depth=N to increase maximum nesting level. 1 error generated.. ', &quot;FunctionGraph(MakeVector{dtype='float64'}(V0, &lt;TensorType(float64, scalar)&gt;, &lt;TensorType(float64, scalar)&gt;,..., &lt;TensorType(float64, scalar)&gt;, &lt;TensorType(float64, scalar)&gt;))&quot;)
  </code></pre>
  <p>It works when the data is 10 timesteps long, so I think the large dataset is causing the problem. Here's my model:</p>
  <pre><code>daily = [insert time series here; length 600]
  mu = pm.Normal(&quot;mu&quot;, mu=0.01, sigma=0.07)
  mu_y = pm.Normal(&quot;mu_y&quot;, mu=0, sigma=1)
  sig_y = pm.HalfNormal(&quot;sig_y&quot;, sigma=1)
  lamb = pm.Beta(&quot;lambda&quot;, alpha=2, beta=40)
  sig_v = pm.TruncatedNormal(&quot;sig_v&quot;, mu=0.5, sigma=0.2, lower=0)
  rho_j = pm.TruncatedNormal(&quot;rho_j&quot;, mu=0, sigma=0.3, lower=-1, upper=1)
  mu_v = pm.HalfNormal(&quot;mu_v&quot;, sigma=1)
  rho = pm.TruncatedNormal(&quot;rho&quot;, mu=0, sigma=0.3, lower=-1, upper=1)
  
  epsSigma = tt.stack([1.0, rho, rho, 1.0]).reshape((2, 2))
  eps = [pm.MvNormal(&quot;eps&quot;+str(i), mu=np.zeros(2), cov=epsSigma, shape=2) for i in range(len(daily))]
  alphabeta = pm.MvNormal(&quot;alphabeta&quot;, mu=np.zeros(2), cov=np.eye(2), shape=2)
  
  Zv = pm.Exponential(&quot;Zv&quot;, lam=mu_v, dims=&quot;date&quot;)
  
  Zy_mu = mu_y + rho_j * Zv
  Zy = pm.Normal('Zy', mu=Zy_mu, sigma=sig_y, dims='date')
  
  J = pm.Bernoulli(&quot;J&quot;, p=lamb, dims=&quot;date&quot;)
  
  
  # this is where the problems start
  
  
  
  V = [i for i in range(len(daily))]
  V[0] = pm.TruncatedNormal(&quot;V0&quot;, mu=alphabeta[0], sigma=1, lower=0)
  for t in range(1, len(V)):
      V[t] = alphabeta[0] + alphabeta[1] * V[t-1]
      V[t] = V[t] + sig_v * eps[t][1]
      V[t] = V[t] * pm.math.sqrt(10**-8 + pm.math.maximum(V[t-1], 0)) 
      V[t] = V[t] + J[t] * Zv[t]
  V = pm.Normal(&quot;V&quot;, mu=tt.stack(V), sigma=0.05, dims=&quot;date&quot;)
  
  Y = [i for i in range(len(daily))]
  Y[0] = pm.Normal(&quot;YO&quot;, mu=mu, sigma=0.1)
  for t in range(1, len(Y)):
      Y[t] = mu + pm.math.sqrt(10**-8 + pm.math.maximum(V[t-1], 0)) * eps[t][1] + J[t] * Zy[t]
  Y_obs = pm.Normal(&quot;Y_obs&quot;, mu=tt.stack(Y), sigma=0.05, dims=&quot;date&quot;, observed=daily_obs)
  </code></pre>
  <p>From <a href="https://stackoverflow.com/questions/53670475/bayesian-linear-regression-with-pymc3-and-a-large-dataset-bracket-nesting-leve">this post</a>, it seems that the issue is with the for loops used to define the length-600 vectors <code>V</code> and <code>Y</code>, but it doesn't seem to be possible to do the definitions without loops, given the nature of the recurrence relations (here they are in LaTeX for readability):</p>
  <p><a href="https://i.stack.imgur.com/TOxZr.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/TOxZr.png" alt="enter image description here" /></a></p>
  <p>Is there a better way to define these variables?</p>