diff --git a/GPy/likelihoods/gaussian.py b/GPy/likelihoods/gaussian.py
index 7b001e17..ef4b26b1 100644
--- a/GPy/likelihoods/gaussian.py
+++ b/GPy/likelihoods/gaussian.py
@@ -326,7 +326,7 @@ class Gaussian(Likelihood):
         dF_dtheta = -0.5/lik_var + 0.5*(np.square(Y) + np.square(m) + v - 2*m*Y)/(lik_var**2)
         return F, dF_dmu, dF_dv, dF_dtheta.reshape(1, Y.shape[0], Y.shape[1])
 
-class Heteroscedastic_Gaussian(Gaussian):
+class HeteroscedasticGaussian(Gaussian):
     def __init__(self, Y_metadata, gp_link=None, variance=1., name='het_Gauss'):
         if gp_link is None:
             gp_link = link_functions.Identity()
@@ -335,7 +335,7 @@ class Heteroscedastic_Gaussian(Gaussian):
             print("Warning, Exact inference is not implemeted for non-identity link functions,\
             if you are not already, ensure Laplace inference_method is used")
 
-        super(Heteroscedastic_Gaussian, self).__init__(gp_link, np.ones(Y_metadata['output_index'].shape[0])*variance, name)
+        super(HeteroscedasticGaussian, self).__init__(gp_link, np.ones(Y_metadata['output_index'].shape[0])*variance, name)
 
     def exact_inference_gradients(self, dL_dKdiag,Y_metadata=None):
         return dL_dKdiag[Y_metadata['output_index']][:,0]
@@ -351,4 +351,4 @@ class Heteroscedastic_Gaussian(Gaussian):
                 var += np.atleast_3d(np.eye(var.shape[0])*self.variance)
         else:
             var += self.variance
-        return mu, var
\ No newline at end of file
+        return mu, var
diff --git a/GPy/models/gp_heteroscedastic_regression.py b/GPy/models/gp_heteroscedastic_regression.py
index 19cb18d8..63c6352a 100644
--- a/GPy/models/gp_heteroscedastic_regression.py
+++ b/GPy/models/gp_heteroscedastic_regression.py
@@ -16,6 +16,8 @@ class GPHeteroscedasticRegression(GP):
     :param X: input observations
     :param Y: observed values
     :param kernel: a GPy kernel, defaults to rbf
+
+    NB: This model does not make inference on the noise outside the training set
     """
     def __init__(self, X, Y, kernel=None, Y_metadata=None):
 
@@ -30,10 +32,7 @@ class GPHeteroscedasticRegression(GP):
             kernel = kern.RBF(X.shape[1])
 
         #Likelihood
-        #likelihoods_list = [likelihoods.Gaussian(name="Gaussian_noise_%s" %j) for j in range(Ny)]
-        noise_terms = np.unique(Y_metadata['output_index'].flatten())
-        likelihoods_list = [likelihoods.Gaussian(name="Gaussian_noise_%s" %j) for j in noise_terms]
-        likelihood = likelihoods.MixedNoise(likelihoods_list=likelihoods_list)
+        likelihood = likelihoods.HeteroscedasticGaussian(Y_metadata)
 
         super(GPHeteroscedasticRegression, self).__init__(X,Y,kernel,likelihood, Y_metadata=Y_metadata)