Problems when fitting a reinforcement learning model to data

Questions
4 
with pm.Model() as RW_model:
    ##############################################
    beta0 = pm.Normal('beta0',0,1)
    beta1 = pm.Normal('beta1',0,1)
    ##############################################
    lr = pm.Beta('lr', 1,1, shape=no_subjects)
    ##############################################
    
    # step function
    vS = 0.5 * tt.ones((no_subjects,2), dtype='float64') 
    combined_v = 0.5 * tt.ones((no_subjects,1), dtype='float64') 

    outputs, updates = theano.scan(
        fn=update_RW,
        sequences=[stim_matrix, shock_matrix],
        outputs_info=[vS, combined_v],
        non_sequences=[lr, no_subjects])
    
    combined_v_ = outputs[1].reshape(est_scr_matrix.shape)
    est_scr_model = combined_v_*beta1 +beta0
    ##############################################
    # Model error
    eps = pm.HalfNormal('eps', 0.5)
    
    # Data likelihood
    scrs = pm.Normal('scrs',est_scr_model,eps,observed=est_scr_matrix)
    ##############################################
    
    trace = pm.sample(target_accept=0.9, chains=4, cores=10, return_inferencedata=True)

The snippet above contains the simplest version of my model. I’ve narrowed down the problem to the beta0 parameter and the equation est_scr_model = combined_v_*beta1 +beta0. The divide by zero warning appears when the following exception is thrown:

".../pymc3/step_methods/hmc/quadpotential.py", line 272, in raise_ok
    raise ValueError("\n".join(errmsg))
ValueError: Mass matrix contains zeros on the diagonal. 
The derivative of RV `beta0`.ravel()[0] is zero.

It looks like something is wrong when quadpotential is performed to calculate beta0. Strangely, the problem is gone if I only keep one of the betas, i.e. making the equation either est_scr_model = combined_v_*beta1 or est_scr_model = combined_v_ + beta0 solves the problem.

What should I do if I want to keep both betas?