I’m running into an issue with Theano when building a custom hierarchical PyMC3 model using a black box likelihood. When running the model with too many subjects, I receive the following error:

'Compilation failed (return status=1): /Users/matthewmurrow/.theano/compiledir_macOS-10.16-x86_64-i386-64bit-i386-3.8.8-64/tmplbcx0v1h/mod.cpp:35298:32: fatal error: bracket nesting level exceeded maximum of 256. if (!PyErr_Occurred())

The relevant part of the code is shown below. N_subjects is the number of individuals with data I am including in the model. The logl_fun callouts in the pm.Potentials reference the custom likelihood function itself. From looking at other posts, the issue seems to be related to the python loop I’m doing in the PyMC function. Is this correct, and if so, will a Theano scan loop solve the problem?

with pm.Model() as hm:

`# prior distributions mean_mu = pm.Normal('mean_mu0', mu = 1.0, sigma = 2.0) sd_mu = pm.Gamma('sd_mu0', alpha = 1.0, beta = 1.0) mean_mu1 = pm.Normal('mean_mu1', mu = 1.0, sigma = 1.0) sd_mu1 = pm.Gamma('sd_mu1', alpha = 1.0, beta = 1.0) tnd = pm.Uniform("tnd", lower = 0.0, upper = 0.5) w = pm.Uniform("w", lower = 0.3, upper = 0.7) b = pm.Uniform("b", lower = 0.1, upper = 1.0) # group level distributions mu = pm.Normal('mu_c0', mu = mean_mu, sigma = sd_mu, shape = N_subjects) mu1 = pm.Normal('mu_c1', mu = mean_mu1, sigma = sd_mu1, shape = N_subjects) psi_temp = [0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] psi_temp[0] = w psi_temp[2] = tnd psi_pymc = tt.as_tensor_variable(psi_temp) for ii in range(N_subjects): phi_pymc[ii] = tt.as_tensor_variable([mu[ii], 0.0, 0.0, 1.0, b, 0.0, 0.0, np.log10(1.0), np.log10(1.0), np.log10(1.0), 1.0, 0.0, 0.0, 1.0]) phi_pymc1[ii] = tt.as_tensor_variable([mu1[ii], 0.0, 0.0, 1.0, b, 0.0, 0.0, np.log10(1.0), np.log10(1.0), np.log10(1.0), 1.0, 0.0, 0.0, 1.0]) pm.Potential('c0' + str(ii), logl_fun(phi_pymc[ii], psi_pymc, ii)) pm.Potential('c1' + str(ii), logl_fun1(phi_pymc1[ii], psi_pymc, ii)) step = pm.DEMetropolisZ(tune = 'lambda') samples = 50000 tune = 0 chains = 2 cores = 2 idata = pm.sample(draws = samples, tune = tune, step = step, compute_convergence_checks = False, chains = chains, cores = cores, return_inferencedata = True, discard_tuned_samples = False) idata.to_netcdf('test' + '.nc')`