# Blockwise Metropolis

**URL:** <https://discourse.pymc.io/t/blockwise-metropolis/8784>\
**Category:** Questions\
**Created:** [February 11, 2022, 4:11pm UTC](https://discourse.pymc.io/t/blockwise-metropolis/8784 "2022-02-11T16:11:41Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![aarcher](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/aarcher/32/4704_2.png) [@aarcher](https://discourse.pymc.io/u/aarcher)\
**Post date:** [February 11, 2022, 4:11pm UTC](https://discourse.pymc.io/t/blockwise-metropolis/8784/1 "2022-02-11T16:11:41Z")

</div>

Hi, I’m new to the pyMC community. I noticed that pm.sample(draws=10000, step=pm.Metropolis()) did a Metropolis within Gibbs sampling across all parameters. Is it possible to have pm.Metropolis() do blockwise sampling using a multivariate normal distribution? I’ve tried doing a compound step with disjoint parameter partitions using pm.Metropolis(). But this resulted in a Metropolis within Gibbs sampling across all parameters.

---

<div class="post-metadata">

**Author:** ![ckrapu](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ckrapu/32/1366_2.png) [@ckrapu](https://discourse.pymc.io/u/ckrapu)\
**Post date:** [February 11, 2022, 4:26pm UTC](https://discourse.pymc.io/t/blockwise-metropolis/8784/2 "2022-02-11T16:26:51Z")

</div>

If you use a `MvNormal` distribution, the default sampler should attempt to update all elements of that random variable simultaneously.

That said, my experience is that using a multivariate Metropolis proposal without NUTS or HMC is usually pretty poor in performance and scales exponentially poorly with the dimension of the RV. The default sampler isn’t going to use blockwise **conjugate** updates for a multivariate normal, which is how I’ve often seen handwritten MCMC implementations do it.

---

<div class="post-metadata">

**Author:** ![aarcher](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/aarcher/32/4704_2.png) [@aarcher](https://discourse.pymc.io/u/aarcher)\
**Post date:** [February 11, 2022, 5:05pm UTC](https://discourse.pymc.io/t/blockwise-metropolis/8784/3 "2022-02-11T17:05:19Z")

</div>

Thank you for your response!

I’m using my own custom likelihood function. My code is as follows:

```auto
class LogLike(at.Op):

    itypes = [at.dvector] # expects a vector of parameter values when called
    otypes = [at.dscalar] # outputs a single scalar value (the log likelihood)

    def __init__ (self, loglik, rtol, atol):
        """
        Initialise the Op with various things that our log-likelihood function
        requires. 

        Parameters
        ----------
        loglik:
            The log-likelihood
        """

        # add inputs as class attributes
        self.likelihood = lambda sample: loglik(sample)

    def perform(self, node, inputs, outputs):
        # the method that is used when calling the Op
        (theta,) = inputs # this will contain my variables
        # call the log-likelihood function
        logl = self.likelihood(theta)

logl = LogLike(loglik)
with pm.Model():
  variables = []
  for param_name in variable_name:
     variables.append(pm.Normal(param_name, mu=0., sigma=1.))
  theta = at.as_tensor_variable(variables)
  # use a Potential to "call" the Op and include it in the logp computation
  pm.Potential("likelihood", logl(theta))
  mh = pm.sample(draws=int(nsamples), step=pm.Metropolis())

```

To clarify, are you recommending I change

```auto
variables = []
for param_name in param_names:
   variables.append(pm.Normal(param_name, mu=0., sigma=1.))

```

to

`variables = MvNormal("variables", mu=np.zeros(len(param_names)),cov(np.diag(ones(len(param_names))))?`

---

<div class="post-metadata">

**Author:** ![ckrapu](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ckrapu/32/1366_2.png) [@ckrapu](https://discourse.pymc.io/u/ckrapu)\
**Post date:** [February 11, 2022, 8:38pm UTC](https://discourse.pymc.io/t/blockwise-metropolis/8784/4 "2022-02-11T20:38:10Z")

</div>

Yes, that’s precisely right. Even though it’s the same underlying probability model, this way will force the sampler logic to treat the normal variates as a block.