Hey ,

I am pretty new to pymc and think I am missing something obvious. I have an experiment to measure the impact of push notifications on completing a task. The treatment is a binary 1 to receive notifications and 0 to not.

I can get the treatment effect with the slope form a simple linear regression:

```
Normal('slope' ....)
Normal('intercept' ....)
Normal('y', slope * data.treatment + intercept, observed=data.y)
```

In reality though about 60% of users have notification turned off, so even if they are assigned to the notification group they do not get any. I was thinking prob programming could be really cool here if I modeled the uncertainty with a Bernoulli distribution like:

```
prob_of_notifications = Normal('prob_of_notifications', .6, .1)
b = Bernoulli('b', p=prob_of_notifications, size=len(data))
real_treatment = Deterministic('real_treatment', pm.math.switch(data.treatment, b, data.treatement))
pm.Normal('y', slope * real_treatment + intercept, observed=data.y)
```

When I do this though its becomes way too slow. I tried replacing with at.random.normal etc for those but would not let me do that. Does this type of model make any sense and is there any way to make it efficient if I don’t care about generating a posterior for the real_treatment?

(also is this an appropriate place to ask this type of question? sorry new here!)