I am trying to implement model checks for a Bayesian beta binomial model using the workflow suggested by Betancourt (https://betanalpha.github.io/assets/case_studies/principled_bayesian_workflow.html#51_first_iteration:_a_journey_of_a_thousand_models_begins_with_a_single_step).

I am currently struggling computing the prior_sd of the parameter **p** (which is required to compute the posterior shrinkage). I know how to compute the sd of a beta distribution based on parameters alpha and beta but I don’t know how to do it when alpha and beta are not numbers but distributions. Could anyone help or point me in the right direction? Below is the current implementation of the model (in pymc3).

with pm.Model() as beta_binomial:

```
# Priors
intercept = pm.Normal('intercept', mu=np.log(0.1), sd=0.05)
beta = pm.Normal('beta', mu=np.log(0.75), sd=0.05)
# Beta parameters
p_mu = pm.Deterministic('p_mu', pm.math.invlogit(intercept + beta * X_predictor))
p_kappa = pm.HalfNormal('p_kappa', 10)
p_alpha = p_mu * p_kappa
p_beta = (1-p_mu) * p_kappa
# Outcome definition
p = pm.Beta('p', alpha=p_alpha, beta=p_beta, shape=X_predictor.shape[0])
Y = pm.Binomial('num_successes', n=num_trials, p=p, observed=num_successes)
```