Added compute of jacobian.

GPy/core/parameterization/transformations.py

@@ -31,6 +31,13 @@ class Transformation(object):
         raise NotImplementedError
     def finv(self, model_param):
         raise NotImplementedError
+    def jacobianfactor(self, model_param):
+        """
+        Return log|det J| where J is the Jacobian of the inverse of the
+        transformation.
+        """
+        return np.abs([self.gradfactor(np.array([theta]), np.ones((1,)))[0]
+                       for theta in model_param])
     def gradfactor(self, model_param, dL_dmodel_param):
         """ df(opt_param)_dopt_param evaluated at self.f(opt_param)=model_param, times the gradient dL_dmodel_param,
 
@@ -67,7 +74,7 @@ class Logexp(Transformation):
         return np.where(x>_lim_val, x, np.log1p(np.exp(np.clip(x, -_lim_val, _lim_val)))) + epsilon
         #raises overflow warning: return np.where(x>_lim_val, x, np.log(1. + np.exp(x)))
     def finv(self, f):
-        return np.where(f>_lim_val, f, np.log(np.exp(f+1e-20) - 1.))
+        return np.where(f>_lim_val, f, np.log(np.exp(f+epsilon) - 1.))
     def gradfactor(self, f, df):
         return np.einsum('i,i->i', df, np.where(f>_lim_val, 1., 1. - np.exp(-f)))
     def initialize(self, f):

ib_tests/test_transformation_of_pdf.py

@@ -0,0 +1,57 @@
+"""
+Test the transformation of a PDF.
+
+Author:
+    Ilias Bilionis
+
+Date:
+    8/4/2015
+
+"""
+
+
+import sys
+import os
+# Make sure we load the GP that is here
+sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..')))
+import GPy
+import matplotlib.pyplot as plt
+import numpy as np
+import scipy.stats as st
+import scipy.integrate as integrate
+
+
+if __name__ == '__main__':
+    p_theta = st.lognorm(.9)
+
+    # Plot the PDF of theta
+    fig, ax = plt.subplots()
+    theta = np.linspace(0.0001, 8., 100)
+    ax.plot(theta, p_theta.pdf(theta), linewidth=2, label='True PDF')
+    ax.set_xlabel(r'$\theta$', fontsize=16)
+    ax.set_ylabel(r'$p(\theta)$', fontsize=16)
+    
+    # Now let's look at the transformation phi = log(exp(theta - 1))
+    t = GPy.constraints.Logexp()
+    t.plot()
+    # Plot the transformed probability density
+    phi = np.linspace(-8, 8, 100)
+    fig, ax = plt.subplots()
+    ax.plot(phi, p_theta.pdf(t.f(phi)) * t.jacobianfactor(t.f(phi)), linewidth=2,
+            label='Transformed PDF')
+    # Now find the normalization constant for the naive transformation of the
+    # PDF
+    p_phi_prop = lambda(phi): p_theta.pdf(t.f(phi))
+    c = integrate.quad(p_phi_prop, -np.inf, np.inf)[0]
+    p_phi = lambda(phi): p_phi_prop(phi) / c
+    ax.plot(phi, p_phi(phi), '--', linewidth=2, label='Naively transformed PDF')
+    # Now let's draw some samples of theta and transform them so that we see
+    # which one is right
+    theta_s = p_theta.rvs(100000)
+    phi_s = t.finv(theta_s)
+    ax.hist(phi_s, normed=True, bins=100, alpha=0.25, label='Empirical')
+    ax.set_xlim(-3, 10)
+    ax.set_xlabel(r'transformed $\theta$', fontsize=16)
+    ax.set_ylabel('PDF', fontsize=16)
+    plt.legend(loc='best')
+    plt.show(block=True)