Sometimes when working with simple models (1-4 parameters), I think it would be cool to have the option to do grid estimation (e.g., to debug / compare with MCMC in case of divergences).

Obviously it is not very difficult to write your own grid estimator, but it would be interesting if we could use the PyMC syntax to set up the model as usual and just call pm.grid(points=dict()) or something similar. You would just need to specify the points in the grid that you are interested in computing. Here is a pseudo-example of how this could look like:

```
with pm.Model() as m:
b0 = pm.Normal('b0', mu=0, sigma=1)
b1 = pm.Normal('b1', mu=0, sigma=1)
sigma = pm.HalfNormal('sigma', sigma=10)
like = pm.Normal('like', mu=b0 + x@b1, sigma=sigma, observed=y)
pm.grid({'b0': np.linspace(-5, 5, 100),
'b1': np.linspace(-5, 5, 100);
'sigma':np.linspace(0.001, 50, 50)})
```

This would make it trivial to do grid sampling for people familiar with the PyMC syntax (and take advantage of the myriad of custom built-in distributions), and one would always be a single line away from being able to run MCMC with the exact same model.

What do you think?