diff --git a/sklearn/gaussian_process/gaussian_process.py b/sklearn/gaussian_process/gaussian_process.py
index 2ab890592ce8e10c60d903f956f3d3c0d990dc4e..5470c809b69d0780a0b9d247837e732e9f222250 100644
--- a/sklearn/gaussian_process/gaussian_process.py
+++ b/sklearn/gaussian_process/gaussian_process.py
@@ -701,9 +701,9 @@ class GaussianProcess(BaseEstimator, RegressorMixin):
 
             constraints = []
             for i in range(self.theta0.size):
-                constraints.append(lambda log10t:
+                constraints.append(lambda log10t, i=i:
                                    log10t[i] - np.log10(self.thetaL[0, i]))
-                constraints.append(lambda log10t:
+                constraints.append(lambda log10t, i=i:
                                    np.log10(self.thetaU[0, i]) - log10t[i])
 
             for k in range(self.random_start):
diff --git a/sklearn/gaussian_process/tests/test_gaussian_process.py b/sklearn/gaussian_process/tests/test_gaussian_process.py
index 61d617690edb978d20d78de6ec298211cb2a7a56..ef31237bb4571d461855f822451bf16892430b66 100644
--- a/sklearn/gaussian_process/tests/test_gaussian_process.py
+++ b/sklearn/gaussian_process/tests/test_gaussian_process.py
@@ -58,15 +58,21 @@ def test_2d(regr=regression.constant, corr=correlation.squared_exponential,
                   [-2.87600388, 6.74310541],
                   [5.21301203, 4.26386883]])
     y = g(X).ravel()
+
+    thetaL = [1e-4] * 2
+    thetaU = [1e-1] * 2
     gp = GaussianProcess(regr=regr, corr=corr, beta0=beta0,
-                         theta0=[1e-2] * 2, thetaL=[1e-4] * 2,
-                         thetaU=[1e-1] * 2,
+                         theta0=[1e-2] * 2, thetaL=thetaL,
+                         thetaU=thetaU,
                          random_start=random_start, verbose=False)
     gp.fit(X, y)
     y_pred, MSE = gp.predict(X, eval_MSE=True)
 
     assert_true(np.allclose(y_pred, y) and np.allclose(MSE, 0.))
 
+    assert_true(np.all(gp.theta_ >= thetaL)) # Lower bounds of hyperparameters
+    assert_true(np.all(gp.theta_ <= thetaU)) # Upper bounds of hyperparameters
+
 
 def test_2d_2d(regr=regression.constant, corr=correlation.squared_exponential,
                random_start=10, beta0=None):