Hi PyMC3 team,

I am currently trying to implement the following model of causal inference that was originally demonstrated in Stan here.

Below is the following code I have run. First, I simulate the dataset:

```
N = 500 # Total sample size
alpha = 1.0 # Intercept in the y-model
tau = 0.25 # Treatment effect
# Assignment mechanism
N_t = 200
W = np.random.choice(2, N) # Randomly assign treatment
# Simulate the outcomes
mu_c = alpha # Mean for control
sd_c = 1.0 # Standard deviation for control
mu_t = alpha + tau # Mean for treatment
sd_t = 1.0 # Standard deviation for control
means = [mu_c, mu_t]
rho = 0.0 # Correlation
cov_mat = np.array([[sd_c**2.0, rho * sd_c * sd_t], [rho * sd_c * sd_t, sd_t**2.0]]) # Covariant matrix
outcomes = np.random.multivariate_normal(means, cov_mat, size=N) # Generate via multivariate normal
Y0 = outcomes[:, 0] # potential outcome if W = 1
Y1 = outcomes[:, 1] # potential outcome if W = 0
tau_unit = Y1 - Y0
# The realization of potential outcomes by the assignment mechanism
Y_obs = Y0 * (1.0 - W) + Y1 * W # Observed outcomes
Y_mis = Y0 * W + Y1 * (1.0 - W) # Missing outcomes
```

And this is the code to create the model in PyMC3:

```
# Now let's try and build this model in PyMC3
with pm.Model() as model:
# Define the priors
alpha = pm.Normal('alpha', mu=0.0, sigma=5.0)
tau = pm.Normal('tau', mu=0.0, sigma=5.0)
sigma_c = pm.Normal('sigma_c', mu=0.0, sigma=5.0)
sigma_t = pm.Normal('sigma_t', mu=0.0, sigma=5.0)
mu = alpha + tau * W
sigma = sigma_t * W + sigma_c * (1.0 - W)
# Define the likelihood
y_obs = pm.Normal('Y_obs', mu=mu, sigma=sigma, observed=Y_obs)
# Sample from the posterior
trace = pm.sample(2000, tune=1000)
```

Which produces the following error:

```
---------------------------------------------------------------------------
SamplingError Traceback (most recent call last)
<ipython-input-414-f57626c5907e> in <module>
13
14 # Sample from the posterior
---> 15 trace = pm.sample(2000, tune=1000)
/usr/local/lib/python3.8/site-packages/pymc3/sampling.py 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, jitter_max_retries, return_inferencedata, idata_kwargs, mp_ctx, pickle_backend, **kwargs)
426 start = deepcopy(start)
427 if start is None:
--> 428 check_start_vals(model.test_point, model)
429 else:
430 if isinstance(start, dict):
/usr/local/lib/python3.8/site-packages/pymc3/util.py in check_start_vals(start, model)
235
236 if not np.all(np.isfinite(initial_eval)):
--> 237 raise SamplingError(
238 "Initial evaluation of model at starting point failed!\n"
239 "Starting values:\n{}\n\n"
SamplingError: Initial evaluation of model at starting point failed!
Starting values:
{'alpha': array(0.), 'tau': array(0.), 'sigma_c': array(0.), 'sigma_t': array(0.)}
Initial evaluation results:
alpha -2.53
tau -2.53
sigma_c -2.53
sigma_t -2.53
Y_obs -inf
Name: Log-probability of test_point, dtype: float64
```

Is there any way I can change the initial starting point so that I don’t encounter this error? Or am I doing something else incorrectly? Any advice on this would be greatly appreciated! Thank you!

Best wishes,

Axel.