diff --git a/gtsam/linear/linear.i b/gtsam/linear/linear.i
index fdf156ff9..f5857b0c5 100644
--- a/gtsam/linear/linear.i
+++ b/gtsam/linear/linear.i
@@ -490,6 +490,8 @@ virtual class GaussianConditional : gtsam::JacobianFactor {
   bool equals(const gtsam::GaussianConditional& cg, double tol) const;
   
   // Standard Interface
+  double evaluate(const gtsam::VectorValues& x) const;
+  double logDensity(const gtsam::VectorValues& x) const;
   gtsam::Key firstFrontalKey() const;
   gtsam::VectorValues solve(const gtsam::VectorValues& parents) const;
   gtsam::JacobianFactor* likelihood(
@@ -543,17 +545,20 @@ virtual class GaussianBayesNet {
   bool equals(const gtsam::GaussianBayesNet& other, double tol) const;
   size_t size() const;
 
-  // Standard interface
   void push_back(gtsam::GaussianConditional* conditional);
   void push_back(const gtsam::GaussianBayesNet& bayesNet);
   gtsam::GaussianConditional* front() const;
   gtsam::GaussianConditional* back() const;
 
+  // Standard interface
+  double evaluate(const gtsam::VectorValues& x) const;
+  double logDensity(const gtsam::VectorValues& x) const;
+
   gtsam::VectorValues optimize() const;
-  gtsam::VectorValues optimize(gtsam::VectorValues given) const;
+  gtsam::VectorValues optimize(const gtsam::VectorValues& given) const;
   gtsam::VectorValues optimizeGradientSearch() const;
   
-  gtsam::VectorValues sample(gtsam::VectorValues given) const;
+  gtsam::VectorValues sample(const gtsam::VectorValues& given) const;
   gtsam::VectorValues sample() const;
   gtsam::VectorValues backSubstitute(const gtsam::VectorValues& gx) const;
   gtsam::VectorValues backSubstituteTranspose(const gtsam::VectorValues& gx) const;
diff --git a/python/gtsam/tests/test_GaussianBayesNet.py b/python/gtsam/tests/test_GaussianBayesNet.py
index e4d396cfe..022de8c3f 100644
--- a/python/gtsam/tests/test_GaussianBayesNet.py
+++ b/python/gtsam/tests/test_GaussianBayesNet.py
@@ -29,8 +29,7 @@ def smallBayesNet():
     """Create a small Bayes Net for testing"""
     bayesNet = GaussianBayesNet()
     I_1x1 = np.eye(1, dtype=float)
-    bayesNet.push_back(GaussianConditional(
-        _x_, [9.0], I_1x1, _y_, I_1x1))
+    bayesNet.push_back(GaussianConditional(_x_, [9.0], I_1x1, _y_, I_1x1))
     bayesNet.push_back(GaussianConditional(_y_, [5.0], I_1x1))
     return bayesNet
 
@@ -41,13 +40,21 @@ class TestGaussianBayesNet(GtsamTestCase):
     def test_matrix(self):
         """Test matrix method"""
         R, d = smallBayesNet().matrix()  # get matrix and RHS
-        R1 = np.array([
-            [1.0, 1.0],
-            [0.0, 1.0]])
+        R1 = np.array([[1.0, 1.0], [0.0, 1.0]])
         d1 = np.array([9.0, 5.0])
         np.testing.assert_equal(R, R1)
         np.testing.assert_equal(d, d1)
 
+    def test_sample(self):
+        """Test sample method"""
+        bayesNet = smallBayesNet()
+        sample = bayesNet.sample()
+        self.assertIsInstance(sample, gtsam.VectorValues)
 
-if __name__ == '__main__':
+        # standard deviation is 1.0 for both, so we set tolerance to 3*sigma
+        mean = bayesNet.optimize()
+        self.gtsamAssertEquals(sample, mean, tol=3.0)
+
+
+if __name__ == "__main__":
     unittest.main()