Merge pull request #3057 from bwengals/maunaloa2

A second CO2 example for GPs

Author: bwengals

docs/source/notebooks/GP-MaunaLoa.ipynb

@@ -0,0 +1,196 @@
+"-------------------------------------------------------------------------------------------"
+"  Atmospheric CO2 record based on ice core data before 1958,                               "
+"  (Ethridge et. al., 1996; MacFarling Meure et al., 2006) and                              "
+"  yearly averages of direct observations from Mauna Loa and the South Pole after and       "
+"  including 1958 (from Scripps CO2 Program).                                               "
+"                                                                                           "
+"  in situ data is based on simple average of Mauna Loa and South Pole from values on Jan 1 "
+"                                                                                           "   
+"  For Mauna Loa Observatory, Hawaii: Latitude 19.5�N Longitude 155.6�W Elevation 3397m     "
+"  For South Pole: Latitude 90.0�S Elevation 2810m                                          "
+"                                                                                           "
+"   Scripps CO2 Source: R. F. Keeling, S. C. Piper, A. F. Bollenbacher and S. J. Walker     "
+"   Scripps CO2 Program ( http://scrippsco2.ucsd.edu )                                      "
+"   Scripps Institution of Oceanography (SIO)                                               "
+"   University of California                                                                "
+"   La Jolla, California USA 92093-0244                                                     "
+"                                                                                           "
+"   Status of Scripps CO2 data and correspondence:                                          "
+"                                                                                           "
+"   These data are subject to revision based on recalibration of standard gases. Questions  "
+"   about the data should be directed to Dr. Ralph Keeling ([email protected]) and          "
+"   and Dr. Stephen Piper ([email protected]), Scripps CO2 Program.                  "
+"                                                                                           "
+"                                                                                           "
+"-------------------------------------------------------------------------------------------"
+"                                                                                           "
+Sample date (year), CO2 (ppm)
+  13.3, 276.750
+  29.5, 277.880
+  56.0, 277.380
+ 104.5, 277.510
+ 136.0, 278.130
+ 168.2, 280.050
+ 202.5, 280.720
+ 227.9, 281.490
+ 274.2, 280.130
+ 302.3, 279.830
+ 329.2, 278.920
+ 364.6, 277.050
+ 428.4, 276.910
+ 461.2, 276.720
+ 499.8, 276.360
+ 536.7, 276.000
+ 572.0, 277.560
+ 595.6, 276.900
+ 632.0, 278.250
+ 667.9, 279.390
+ 698.4, 279.700
+ 729.7, 278.530
+ 764.5, 278.530
+ 799.2, 278.550
+ 857.3, 279.340
+ 897.4, 278.910
+ 944.2, 279.120
+ 968.2, 278.460
+1005.0, 280.500
+1025.2, 280.830
+1058.0, 282.760
+1105.4, 282.750
+1159.6, 283.880
+1207.4, 283.600
+1257.6, 282.110
+1275.8, 281.130
+1306.5, 281.500
+1349.7, 280.060
+1411.3, 279.620
+1429.3, 279.540
+1431.0, 282.520
+1560.4, 281.750
+1588.3, 281.030
+1610.4, 271.830
+1628.9, 274.500
+1640.1, 276.620
+1689.6, 276.250
+1722.8, 276.940
+1734.1, 278.230
+1742.7, 276.740
+1752.0, 276.390
+1763.5, 276.320
+1773.7, 277.780
+1780.6, 276.780
+1794.4, 281.540
+1799.6, 281.150
+1814.2, 284.340
+1826.2, 281.280
+1834.5, 283.730
+1838.0, 284.070
+1844.0, 286.500
+1846.0, 284.130
+1849.0, 287.730
+1852.3, 288.570
+1854.0, 287.030
+1859.0, 286.480
+1864.0, 285.410
+1867.0, 285.220
+1869.1, 287.690
+1873.0, 287.170
+1874.0, 290.520
+1884.0, 289.810
+1884.4, 289.010
+1886.0, 290.620
+1889.0, 291.870
+1894.0, 293.840
+1896.0, 298.160
+1899.0, 296.100
+1902.0, 295.320
+1904.0, 295.120
+1909.0, 300.450
+1911.5, 298.390
+1914.0, 300.350
+1918.6, 303.270
+1919.0, 303.550
+1923.0, 303.200
+1923.6, 305.200
+1928.8, 307.770
+1929.0, 305.720
+1933.0, 307.210
+1934.5, 307.820
+1936.5, 308.990
+1938.0, 309.560
+1939.0, 310.940
+1940.0, 311.900
+1941.0, 310.700
+1941.5, 310.300
+1942.0, 311.290
+1943.0, 310.770
+1944.0, 311.590
+1945.0, 309.750
+1946.0, 311.460
+1947.0, 310.750
+1948.0, 310.470
+1949.0, 311.150
+1950.0, 312.550
+1953.0, 312.060
+1954.0, 311.680
+1955.0, 313.640
+1957.0, 314.040
+1958.0, 314.600
+1959.0, 315.370
+1960.0, 316.390
+1961.0, 317.490
+1962.0, 318.070
+1963.0, 318.470
+1964.0, 318.880
+1965.0, 319.290
+1966.0, 320.450
+1967.0, 321.430
+1968.0, 322.060
+1969.0, 323.150
+1970.0, 324.480
+1971.0, 325.400
+1972.0, 326.090
+1973.0, 327.680
+1974.0, 329.190
+1975.0, 329.760
+1976.0, 330.890
+1977.0, 332.040
+1978.0, 333.880
+1979.0, 335.200
+1980.0, 336.900
+1981.0, 338.580
+1982.0, 339.660
+1983.0, 340.890
+1984.0, 342.890
+1985.0, 344.140
+1986.0, 345.540
+1987.0, 346.950
+1988.0, 349.190
+1989.0, 351.060
+1990.0, 352.280
+1991.0, 353.780
+1992.0, 354.860
+1993.0, 355.690
+1994.0, 356.770
+1995.0, 358.680
+1996.0, 360.580
+1997.0, 361.830
+1998.0, 363.630
+1999.0, 366.360
+2000.0, 367.610
+2001.0, 368.910
+2002.0, 370.690
+2003.0, 373.100
+2004.0, 375.190
+2005.0, 376.940
+2006.0, 379.290
+2007.0, 381.010
+2008.0, 383.080
+2009.0, 384.820
+2010.0, 386.720
+2011.0, 388.930
+2012.0, 390.710
+2013.0, 393.310
+2014.0, 395.880
+2015.0, 397.840
+2016.0, 400.660

docs/source/notebooks/GP-MaunaLoa2.ipynb

@@ -16,6 +16,7 @@
            'Periodic',
            'WarpedInput',
            'Gibbs',
+           'ScaledCov',
            'Coregion']
 
 
@@ -542,15 +543,49 @@ def full(self, X, Xs=None):
     def diag(self, X):
         return tt.alloc(1.0, X.shape[0])
 
-def handle_args(func, args):
-    def f(x, args):
-        if args is None:
-            return func(x)
+
+class ScaledCov(Covariance):
+    R"""
+    Construct a kernel by multiplying a base kernel with a scaling
+    function defined using Theano.  The scaling function is
+    non-negative, and can be parameterized.
+
+    .. math::
+       k(x, x') = \phi(x) k_{\text{base}}(x, x') \phi(x')
+
+    Parameters
+    ----------
+    cov_func: Covariance
+        Base kernel or covariance function
+    scaling_func : callable
+        Theano function of X and additional optional arguments.
+    args : optional, tuple or list of scalars or PyMC3 variables
+        Additional inputs (besides X or Xs) to lengthscale_func.
+    """
+    def __init__(self, input_dim, cov_func, scaling_func, args=None, active_dims=None):
+        super(ScaledCov, self).__init__(input_dim, active_dims)
+        if not callable(scaling_func):
+            raise TypeError("scaling_func must be callable")
+        if not isinstance(cov_func, Covariance):
+            raise TypeError("Must be or inherit from the Covariance class")
+        self.cov_func = cov_func
+        self.scaling_func = handle_args(scaling_func, args)
+        self.args = args
+
+    def diag(self, X):
+        X, _ = self._slice(X, None)
+        cov_diag = self.cov_func(X, diag=True)
+        scf_diag = tt.square(tt.flatten(self.scaling_func(X, self.args)))
+        return cov_diag * scf_diag
+
+    def full(self, X, Xs=None):
+        X, Xs = self._slice(X, Xs)
+        scf_x = self.scaling_func(X, self.args)
+        if Xs is None:
+            return tt.outer(scf_x, scf_x) * self.cov_func(X)
         else:
-            if not isinstance(args, tuple):
-                args = (args,)
-            return func(x, *args)
-    return f
+            scf_xs = self.scaling_func(Xs, self.args)
+            return tt.outer(scf_x, scf_xs) * self.cov_func(X, Xs)
 
 
 class Coregion(Covariance):
@@ -615,3 +650,16 @@ def diag(self, X):
         X, _ = self._slice(X, None)
         index = tt.cast(X, 'int32')
         return tt.diag(self.B)[index.ravel()]
+
+
+def handle_args(func, args):
+    def f(x, args):
+        if args is None:
+            return func(x)
+        else:
+            if not isinstance(args, tuple):
+                args = (args,)
+            return func(x, *args)
+    return f
+
+

docs/source/notebooks/GP-MeansAndCovs.ipynb

@@ -552,6 +552,30 @@ def test_raises(self):
             pm.gp.cov.Gibbs(3, lambda x: x, active_dims=[0,1])
 
 
+class TestScaledCov(object):
+    def test_1d(self):
+        X = np.linspace(0, 1, 10)[:, None]
+        def scaling_func(x, a, b):
+            return a + b*x
+        with pm.Model() as model:
+            cov_m52 = pm.gp.cov.Matern52(1, 0.2)
+            cov = pm.gp.cov.ScaledCov(1, scaling_func=scaling_func, args=(2, -1), cov_func=cov_m52)
+        K = theano.function([], cov(X))()
+        npt.assert_allclose(K[0, 1], 3.00686, atol=1e-3)
+        K = theano.function([], cov(X, X))()
+        npt.assert_allclose(K[0, 1], 3.00686, atol=1e-3)
+        # check diagonal
+        Kd = theano.function([], cov(X, diag=True))()
+        npt.assert_allclose(np.diag(K), Kd, atol=1e-5)
+
+    def test_raises(self):
+        cov_m52 = pm.gp.cov.Matern52(1, 0.2)
+        with pytest.raises(TypeError):
+            pm.gp.cov.ScaledCov(1, cov_m52, "str is not callable")
+        with pytest.raises(TypeError):
+            pm.gp.cov.ScaledCov(1, "str is not Covariance object", lambda x: x)
+
+
 class TestHandleArgs(object):
     def test_handleargs(self):
         def func_noargs(x):