Laplace approximation not handling non-scalar parameters

[dar326]

I have been trying to utilize the new laplace function for estimating the posterior distribution using a quadratic approximation. It seems that there is an issue when trying to use non-scalar parameters. This is a slight modification of the example that is provided in the documentation for the function:

y = np.array([2642, 3503, 4358]*10)
cat = np.array([0]*15 + [1]*15)
with pm.Model() as m:
    logsigma = pm.Uniform("logsigma", 1, 100)
    mu = pm.Uniform("mu", -10000, 10000, shape=2)
    yobs = pm.Normal("y", mu=mu[cat], sigma=pm.math.exp(logsigma), observed=y)
    idata = laplace([mu, logsigma], model=m)

The following error is generated:

TypeError: ('Wrong number of dimensions: expected 1, got 0 with shape ().', 'Container name "mu"')

Interestingly, this error occurs even when the shape is not greater than 1 and no indexing is performed:

y = np.array([2642, 3503, 4358]*10)
with pm.Model() as m:
    logsigma = pm.Uniform("logsigma", 1, 100)
    mu = pm.Uniform("mu", -10000, 10000, shape=1)
    yobs = pm.Normal("y", mu=mu, sigma=pm.math.exp(logsigma), observed=y)
    idata = laplace([mu, logsigma], model=m)

This seems to indicate that there is an issue with setting the shape parameter to any value.