diff --git a/GPy/core/parameterization/priors.py b/GPy/core/parameterization/priors.py
index 906ad65f..6269436b 100644
--- a/GPy/core/parameterization/priors.py
+++ b/GPy/core/parameterization/priors.py
@@ -200,37 +200,50 @@ class MultivariateGaussian(Prior):
 
     def __new__(cls, mu=0, var=1):  # Singleton:
         if cls._instances:
-            cls._instances[:] = [instance for instance in cls._instances if instance()]
+            cls._instances[:] = [instance for instance in cls._instances if
+                                 instance()]
             for instance in cls._instances:
-                if np.all(instance().mu == mu) and np.all(instance().var == var):
+                if np.all(instance().mu == mu) and np.all(
+                        instance().var == var):
                     return instance()
-        o = super(Prior, cls).__new__(cls, mu, var)
+        newfunc = super(Prior, cls).__new__
+        if newfunc is object.__new__:
+            o = newfunc(cls)
+        else:
+            o = newfunc(cls, mu, var)
         cls._instances.append(weakref.ref(o))
         return cls._instances[-1]()
 
     def __init__(self, mu, var):
         self.mu = np.array(mu).flatten()
         self.var = np.array(var)
-        assert len(self.var.shape) == 2
-        assert self.var.shape[0] == self.var.shape[1]
+        assert len(self.var.shape) == 2, 'Covariance must be a matrix'
+        assert self.var.shape[0] == self.var.shape[1], \
+            'Covariance must be a square matrix'
         assert self.var.shape[0] == self.mu.size
         self.input_dim = self.mu.size
-        self.inv, self.hld = pdinv(self.var)
-        self.constant = -0.5 * self.input_dim * np.log(2 * np.pi) - self.hld
+        self.inv, _, self.hld, _ = pdinv(self.var)
+        self.constant = -0.5 * (self.input_dim * np.log(2 * np.pi) + self.hld)
+
+    def __str__(self):
+        return 'MultiN(' + str(self.mu) + ', ' + str(np.diag(self.var)) + ')'
 
     def summary(self):
         raise NotImplementedError
 
     def pdf(self, x):
+        x = np.array(x).flatten()
         return np.exp(self.lnpdf(x))
 
     def lnpdf(self, x):
+        x = np.array(x).flatten()
         d = x - self.mu
-        return self.constant - 0.5 * np.sum(d * np.dot(d, self.inv), 1)
+        return self.constant - 0.5 * np.dot(d.T, np.dot(self.inv, d))
 
     def lnpdf_grad(self, x):
+        x = np.array(x).flatten()
         d = x - self.mu
-        return -np.dot(self.inv, d)
+        return - np.dot(self.inv, d)
 
     def rvs(self, n):
         return np.random.multivariate_normal(self.mu, self.var, n)
@@ -247,14 +260,15 @@ class MultivariateGaussian(Prior):
         return self.mu, self.var
 
     def __setstate__(self, state):
-        self.mu = state[0]
+        self.mu = np.array(state[0]).flatten()
         self.var = state[1]
-        assert len(self.var.shape) == 2
-        assert self.var.shape[0] == self.var.shape[1]
+        assert len(self.var.shape) == 2, 'Covariance must be a matrix'
+        assert self.var.shape[0] == self.var.shape[1], \
+            'Covariance must be a square matrix'
         assert self.var.shape[0] == self.mu.size
         self.input_dim = self.mu.size
-        self.inv, self.hld = pdinv(self.var)
-        self.constant = -0.5 * self.input_dim * np.log(2 * np.pi) - self.hld
+        self.inv, _, self.hld, _ = pdinv(self.var)
+        self.constant = -0.5 * (self.input_dim * np.log(2 * np.pi) + self.hld)
 
 def gamma_from_EV(E, V):
     warnings.warn("use Gamma.from_EV to create Gamma Prior", FutureWarning)
diff --git a/GPy/testing/model_tests.py b/GPy/testing/model_tests.py
index 9ee16f99..a241fc9c 100644
--- a/GPy/testing/model_tests.py
+++ b/GPy/testing/model_tests.py
@@ -936,6 +936,34 @@ class GradientTests(np.testing.TestCase):
         #import ipdb;ipdb.set_trace()
         #m.constrain_fixed('.*rbf_var', 1.)
         self.assertTrue(m.checkgrad())
+    
+    def test_simple_MultivariateGaussian_prior(self):
+        X = np.random.multivariate_normal(
+            [1, 5], np.diag([0.5, 0.3]), (100, 1)).reshape(100, 2)
+        Y = X + np.random.randn(100, 2) * 0.05
+        kernel = GPy.kern.RBF(input_dim=2, variance=1,lengthscale=1, ARD=True)
+        kernel.unconstrain()
+        kernel.variance.set_prior(GPy.priors.Gaussian(150, 5))
+        kernel.lengthscale.set_prior(GPy.priors.MultivariateGaussian(
+            np.array([20, 20]), np.diag([5, 5])))
+        m = GPy.models.GPRegression(X, Y, kernel=kernel)
+        m.optimize()
+        print(m.kern.variance)
+        print(m.kern.lengthscale)
+
+    def test_simple_MultivariateGaussian_prior_matrixmean(self):
+        X = np.random.multivariate_normal(
+            [1, 5], np.diag([0.5, 0.3]), (100, 1)).reshape(100, 2)
+        Y = X + np.random.randn(100, 2) * 0.05
+        kernel = GPy.kern.RBF(input_dim=2, variance=1,lengthscale=1, ARD=True)
+        kernel.unconstrain()
+        kernel.variance.set_prior(GPy.priors.Gaussian(150, 5))
+        kernel.lengthscale.set_prior(GPy.priors.MultivariateGaussian(
+            np.array([[20, 20]]), np.diag([5, 5])))
+        m = GPy.models.GPRegression(X, Y, kernel=kernel)
+        m.optimize()
+        print(m.kern.variance)
+        print(m.kern.lengthscale)
 
     def test_multioutput_sparse_regression_1D(self):
         X1 = np.random.rand(500, 1) * 8