Maximum Recursion depth exceeded

My PyMC3 model fails to run pm.sample(...) when K = 10, but for smaller values of K (like 3), it runs fine. Does this throw an error because my model is too big for PyMC3 to handle?
As you can see, I am trying to do Bayesian time series analysis. The model is Markov chain but is also bayesian. This model is a basic one and I intend to do this for very large amount of data
Question 1. Is there anyway so I can avoid the error Max. recursion depth exceeded? and get this running?
I previously got the same error a different but similar model: where i was doing theta[i] = pm.Normal(...) inside an ‘i’ loop (also for other vectorised priors lambda0, lambda1 etc…) which I guess is inefficient. Later, I got rid of the error by vectorising theta and other priors using “shape” and this got rid of the error in my previous model.
Question 2. Is there any way I can vectorise the distributions in my code even more?
Question 3. Why is this error being caused? Is it because its too large for PyMC3 to handle?
Question 4. Is there any other better tool (Python or otherwise but preferably in Python) that is faster for Bayesian time series analysis? I read somewhere that stan/PyStan is faster than PyMC3 at least for Bayesian time series analysis. Is this true? Is there anything faster than stan/PyStan? How does edwardlib compare to PyMC3 and PyStan for the same purpose?

By the way: I need to do this on really large amount of data: ~ 200000 observations would be okay. Running on 10M records would be wonderful: But, I do not know if PyMC (or any other package for that matter) can fit on such a large dataset even if one uses Variational Inference instead of MCMC. Can someone please confirm this?

trace = None
summary_df = None
I = 1
K = 10

with pm.Model() as hotDINA:    
    # Priors: bk, ak, learn_k, ones, ss_k, g_k
    theta   = pm.Normal('theta', mu=0.0, sd=1.0, shape=(I, 1))
    lambda0 = pm.Normal('lambda0', mu=0.0, sd=1.0, shape=(K, 1))    #bk
    lambda1 = pm.Uniform('lambda1', 0.0, 2.5, shape=(K, 1))    #ak
    learn   = pm.Beta('learn', alpha=1, beta=1, shape=(K, 1))
    ones    = pm.Bernoulli('known', p=1, shape=(K, 1))
    ss      = pm.Uniform('ss', 0.5, 1.0, shape=(K, 1))
    g       = pm.Uniform('g', 0, 0.5, shape=(K, 1))
    for i in range(I):
        # t = 0
        for k in range(K):
            prob[i][0][k] = pm.math.invlogit((1.7) * lambda1[k,0] * (theta[i,0] - lambda0[k,0]))
            alpha_name = 'alpha[' + str(i) + ',0,' + str(k) + ']'
            alpha[i][0][k] = pm.Bernoulli(alpha_name, prob[i][0][k])
        for s in range(MAXSKILLS):
            idx = int(idxY[i][0][s] - 1)
            if idx >= K: continue
            py[i][0][idx] = pow(ss[idx,0], alpha[i][0][idx]) * pow(g[idx,0], (1-alpha[i][0][idx]))
        # t = 1,2...T[i]-1
        for t in tqdm(range(1, T[i])):
            for k in range(K):
                alpha[i][t][k] = pm.math.switch(alpha[i][t-1][k], ones[k,0], learn[k,0])
            for s in range(MAXSKILLS):
                idx = int(idxY[i][t][s] - 1)
                if idx >= K: continue
                py[i][t][idx] = pow(ss[idx,0], alpha[i][t][idx]) * pow(g[idx,0], (1-alpha[i][t][idx]))
        start = time.time()
        for t in tqdm(range(T[i])):
            for k in range(K):
                if py[i][t][k] != 0:
                    Y[i][t][k] = pm.Bernoulli(f'y_{i}_{t}_{k}', p=py[i][t][k], observed=[1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0])
        end = time.time()
        res = (end - start)/T[i]
        print("TIME w/o mini batch: ", res)

    trace = pm.sample(1000, tune=500)
    pm.save_trace(trace=trace, directory=".pymc_1.trace", overwrite=True)
    pm.traceplot(trace, ['theta', 'lambda0', 'lambda1'])
    summary_df = pm.stats.summary(trace)


RecursionError                            Traceback (most recent call last)
<ipython-input-36-c58d2f4da611> in <module>
     49     print("DONE")
---> 50     trace = pm.sample(1000, tune=500)
     51     pm.save_trace(trace=trace, directory=".pymc_1.trace", overwrite=True)
     52     pm.traceplot(trace, ['theta', 'lambda0', 'lambda1'])

D:\anaconda3\lib\site-packages\pymc3\ in sample(draws, step, init, n_init, start, trace, chain_idx, chains, cores, tune, progressbar, model, random_seed, discard_tuned_samples, compute_convergence_checks, callback, return_inferencedata, idata_kwargs, **kwargs)
    520         _print_step_hierarchy(step)
    521         try:
--> 522             trace = _mp_sample(**sample_args)
    523         except pickle.PickleError:
    524             _log.warning("Could not pickle model, sampling singlethreaded.")

D:\anaconda3\lib\site-packages\pymc3\ in _mp_sample(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, callback, **kwargs)
   1411     sampler = ps.ParallelSampler(
-> 1412         draws, tune, chains, cores, random_seed, start, step, chain, progressbar
   1413     )
   1414     try:

D:\anaconda3\lib\site-packages\pymc3\ in __init__(self, draws, tune, chains, cores, seeds, start_points, step_method, start_chain_num, progressbar)
    372                 draws, tune, step_method, chain + start_chain_num, seed, start
    373             )
--> 374             for chain, seed, start in zip(range(chains), seeds, start_points)
    375         ]

D:\anaconda3\lib\site-packages\pymc3\ in <listcomp>(.0)
    372                 draws, tune, step_method, chain + start_chain_num, seed, start
    373             )
--> 374             for chain, seed, start in zip(range(chains), seeds, start_points)
    375         ]

D:\anaconda3\lib\site-packages\pymc3\ in __init__(self, draws, tune, step_method, chain, seed, start)
    257         )
    258         try:
--> 259             self._process.start()
    260         except IOError as e:
    261             # Something may have gone wrong during the fork / spawn

D:\anaconda3\lib\multiprocessing\ in start(self)
    110                'daemonic processes are not allowed to have children'
    111         _cleanup()
--> 112         self._popen = self._Popen(self)
    113         self._sentinel = self._popen.sentinel
    114         # Avoid a refcycle if the target function holds an indirect

D:\anaconda3\lib\multiprocessing\ in _Popen(process_obj)
    221     @staticmethod
    222     def _Popen(process_obj):
--> 223         return _default_context.get_context().Process._Popen(process_obj)
    225 class DefaultContext(BaseContext):

D:\anaconda3\lib\multiprocessing\ in _Popen(process_obj)
    320         def _Popen(process_obj):
    321             from .popen_spawn_win32 import Popen
--> 322             return Popen(process_obj)
    324     class SpawnContext(BaseContext):

D:\anaconda3\lib\multiprocessing\ in __init__(self, process_obj)
     87             try:
     88                 reduction.dump(prep_data, to_child)
---> 89                 reduction.dump(process_obj, to_child)
     90             finally:
     91                 set_spawning_popen(None)

D:\anaconda3\lib\multiprocessing\ in dump(obj, file, protocol)
     58 def dump(obj, file, protocol=None):
     59     '''Replacement for pickle.dump() using ForkingPickler.'''
---> 60     ForkingPickler(file, protocol).dump(obj)
     62 #

RecursionError: maximum recursion depth exceeded

Thanks in advance! :slight_smile: