9
9
from nose .tools import assert_true
10
10
11
11
import numpy as np
12
+ import scipy .stats as ss
13
+
14
+ import pylab as pl
12
15
13
16
from sklearn .gaussian_process import GaussianProcess
14
17
from sklearn .gaussian_process import regression_models as regression
@@ -39,7 +42,7 @@ def test_1d(regr=regression.constant, corr=correlation.squared_exponential,
39
42
and np .allclose (MSE2 , 0. , atol = 10 ))
40
43
41
44
42
- def test_2d ( regr = regression .constant , corr = correlation .squared_exponential ,
45
+ def check2d ( X , regr = regression .constant , corr = correlation .squared_exponential ,
43
46
random_start = 10 , beta0 = None ):
44
47
"""
45
48
MLE estimation of a two-dimensional Gaussian Process model accounting for
@@ -49,14 +52,7 @@ def test_2d(regr=regression.constant, corr=correlation.squared_exponential,
49
52
"""
50
53
b , kappa , e = 5. , .5 , .1
51
54
g = lambda x : b - x [:, 1 ] - kappa * (x [:, 0 ] - e ) ** 2.
52
- X = np .array ([[- 4.61611719 , - 6.00099547 ],
53
- [4.10469096 , 5.32782448 ],
54
- [0.00000000 , - 0.50000000 ],
55
- [- 6.17289014 , - 4.6984743 ],
56
- [1.3109306 , - 6.93271427 ],
57
- [- 5.03823144 , 3.10584743 ],
58
- [- 2.87600388 , 6.74310541 ],
59
- [5.21301203 , 4.26386883 ]])
55
+
60
56
y = g (X ).ravel ()
61
57
gp = GaussianProcess (regr = regr , corr = corr , beta0 = beta0 ,
62
58
theta0 = [1e-2 ] * 2 , thetaL = [1e-4 ] * 2 ,
@@ -68,7 +64,7 @@ def test_2d(regr=regression.constant, corr=correlation.squared_exponential,
68
64
assert_true (np .allclose (y_pred , y ) and np .allclose (MSE , 0. ))
69
65
70
66
71
- def test_2d_2d ( regr = regression .constant , corr = correlation .squared_exponential ,
67
+ def check2d_2d ( X , regr = regression .constant , corr = correlation .squared_exponential ,
72
68
random_start = 10 , beta0 = None ):
73
69
"""
74
70
MLE estimation of a two-dimensional Gaussian Process model accounting for
@@ -79,14 +75,7 @@ def test_2d_2d(regr=regression.constant, corr=correlation.squared_exponential,
79
75
b , kappa , e = 5. , .5 , .1
80
76
g = lambda x : b - x [:, 1 ] - kappa * (x [:, 0 ] - e ) ** 2.
81
77
f = lambda x : np .vstack ((g (x ), g (x ))).T
82
- X = np .array ([[- 4.61611719 , - 6.00099547 ],
83
- [4.10469096 , 5.32782448 ],
84
- [
F438
0.00000000 , - 0.50000000 ],
85
- [- 6.17289014 , - 4.6984743 ],
86
- [1.3109306 , - 6.93271427 ],
87
- [- 5.03823144 , 3.10584743 ],
88
- [- 2.87600388 , 6.74310541 ],
89
- [5.21301203 , 4.26386883 ]])
78
+
90
79
y = f (X )
91
80
gp = GaussianProcess (regr = regr , corr = corr , beta0 = beta0 ,
92
81
theta0 = [1e-2 ] * 2 , thetaL = [1e-4 ] * 2 ,
@@ -112,26 +101,128 @@ def test_more_builtin_correlation_models(random_start=1):
112
101
all_corr = ['absolute_exponential' , 'squared_exponential' , 'cubic' ,
113
102
'linear' ]
114
103
104
+ X = np .array ([[- 4.61611719 , - 6.00099547 ],
105
+ [4.10469096 , 5.32782448 ],
106
+ [0.00000000 , - 0.50000000 ],
107
+ [- 6.17289014 , - 4.6984743 ],
108
+ [1.3109306 , - 6.93271427 ],
109
+ [- 5.03823144 , 3.10584743 ],
110
+ [- 2.87600388 , 6.74310541 ],
111
+ [5.21301203 , 4.26386883 ]])
112
+
115
113
for corr in all_corr :
116
114
test_1d (regr = 'constant' , corr = corr , random_start = random_start )
117
- test_2d ( regr = 'constant' , corr = corr , random_start = random_start )
118
- test_2d_2d ( regr = 'constant' , corr = corr , random_start = random_start )
115
+ check2d ( X , regr = 'constant' , corr = corr , random_start = random_start )
116
+ check2d_2d ( X , regr = 'constant' , corr = corr , random_start = random_start )
119
117
120
118
121
119
def test_ordinary_kriging ():
122
120
"""
123
121
Repeat test_1d and test_2d with given regression weights (beta0) for
124
122
different regression models (Ordinary Kriging).
125
123
"""
124
+
125
+ X = np .array ([[- 4.61611719 , - 6.00099547 ],
126
+ [4.10469096 , 5.32782448 ],
127
+ [0.00000000 , - 0.50000000 ],
128
+ [- 6.17289014 , - 4.6984743 ],
129
+ [1.3109306 , - 6.93271427 ],
130
+ [- 5.03823144 , 3.10584743 ],
131
+ [- 2.87600388 , 6.74310541 ],
132
+ [5.21301203 , 4.26386883 ]])
133
+
126
134
test_1d (regr = 'linear' , beta0 = [0. , 0.5 ])
127
135
test_1d (regr = 'quadratic' , beta0 = [0. , 0.5 , 0.5 ])
128
- test_2d ( regr = 'linear' , beta0 = [0. , 0.5 , 0.5 ])
129
- test_2d ( regr = 'quadratic' , beta0 = [0. , 0.5 , 0.5 , 0.5 , 0.5 , 0.5 ])
130
- test_2d_2d ( regr = 'linear' , beta0 = [0. , 0.5 , 0.5 ])
131
- test_2d_2d ( regr = 'quadratic' , beta0 = [0. , 0.5 , 0.5 , 0.5 , 0.5 , 0.5 ])
136
+ check2d ( X , regr = 'linear' , beta0 = [0. , 0.5 , 0.5 ])
137
+ check2d ( X , regr = 'quadratic' , beta0 = [0. , 0.5 , 0.5 , 0.5 , 0.5 , 0.5 ])
138
+ check2d_2d ( X , regr = 'linear' , beta0 = [0. , 0.5 , 0.5 ])
139
+ check2d_2d ( X , regr = 'quadratic' , beta0 = [0. , 0.5 , 0.5 , 0.5 , 0.5 , 0.5 ])
132
140
133
141
134
142
def test_no_normalize ():
135
143
gp = GaussianProcess (normalize = False ).fit (X , y )
136
144
y_pred = gp .predict (X )
137
145
assert_true (np .allclose (y_pred , y ))
146
+
147
+
148
+ @raises (Exception )
149
+ def test_no_multiple_feature ():
150
+ '''
151
+ check that multiple features are not allowed for non-noisy data
152
+ '''
153
+ X = np .array ([[- 4.61611719 , - 6.00099547 ],
154
+ [4.10469096 , 5.32782448 ],
155
+ [0.00000000 , - 0.50000000 ],
156
+ [- 6.17289014 , - 4.6984743 ],
157
+ [- 6.17289014 , - 4.6984743 ],
158
+ [1.3109306 , - 6.93271427 ],
159
+ [- 5.03823144 , 3.10584743 ],
160
+ [- 2.87600388 , 6.74310541 ],
161
+ [5.21301203 , 4.26386883 ]])
162
+
163
+ check2d (X )
164
+
165
+
166
+ def check2dNoisy (X , regr = regression .constant , corr = correlation .squared_exponential ,
167
+ random_start = 10 , beta0 = None ):
168
+ """
169
+ MLE estimation of a two-dimensional Gaussian Process model accounting for
170
+ anisotropy. Check random start optimization.
171
+
172
+ Test the interpolating property.
173
+ """
174
+ b , kappa , e = 5. , .5 , .1
175
+ g = lambda x : b - x [:, 1 ] - kappa * (x [:, 0 ] - e ) ** 2.
176
+
177
+ y = g (X ).ravel () + np .random .normal (0 , 0.1 , X .shape [0 ])
178
+ gp = GaussianProcess (regr = regr , corr = corr , beta0 = beta0 ,
179
+ theta0 = [1e-2 ] * 2 , thetaL = [1e-4 ] * 2 ,
180
+ thetaU = [1. ] * 2 ,
181
+ random_start = random_start , verbose = False , nugget = 0.1 )
182
+ gp .fit (X , y )
183
+ y_pred , MSE = gp .predict (X , eval_MSE = True )
184
+
185
+ fig = pl .figure ()
186
+ pl .plot (X , f (X ), 'r:' , label = u'$f(x) = x\,\sin(x)$' )
187
+ pl .plot (X , y , 'r.' , markersize = 10 , label = u'Observations' )
188
+ pl .plot (X , y_pred , 'b-' , label = u'Prediction' )
189
+ pl .fill (np .concatenate ([X , X [::- 1 ]]),
190
+ np .concatenate ([y_pred - 1.9600 * MSE ,
191
+ (y_pred + 1.9600 * MSE )[::- 1 ]]),
192
+ alpha = .5 , fc = 'b' , ec = 'None' , label = '95% confidence interval' )
193
+ pl .xlabel ('$x$' )
194
+ pl .ylabel ('$f(x)$' )
195
+ pl .ylim (- 10 , 10 )
196
+ pl .legend (loc = 'upper left' )
197
+ pl .show ()
198
+
199
+ print ((np .abs (y_pred - y ) <= ss .norm .ppf (0.999 , y_pred , np .sqrt (MSE ))))
200
+
201
+ assert_true ((np .abs (y_pred - y ) <= ss .norm .ppf (0.999 , y_pred , np .sqrt (
202
+ MSE ))).all ()) #check that difference between prediction and truth is smaller than 99.9% conf int.
203
+
204
+
205
+ def test_1d_noisy (regr = regression .constant , corr = correlation .absolute_exponential ,
206
+ random_start= 10 , beta0 = None ):
207
+ """
208
+ MLE estimation of a one-dimensional Gaussian Process model.
209
+ Check random start optimization with noisy / duplicate inputs.
210
+
211
+ Test the interpolating property.
212
+ """
213
+
214
+ X = np .atleast_2d ([1. , 3. , 5. , 6. , 7. , 8. , 9. , 10. ] * 2 ).T
215
+ x = np .atleast_2d (np .linspace (0 , 10 , 50 )).T
216
+
217
+ y = f (X ).ravel () + np .random .normal (0 , 0.1 , len (X ))
218
+
219
+ gp = GaussianProcess (regr = regr , corr = corr , beta0 = beta0 ,
220
+ theta0 = 1e-2 , thetaL = 1e-4 , thetaU = 1e-1 ,
221
+ random_start = random_start , verbose = False , nugget = 0.01 ).fit (X , y )
222
+ y_pred , MSE = gp .predict (x , eval_MSE = True )
223
+
224
+ y = f (x ).ravel ()
225
+ assert_true ((np .abs (y_pred - y ) <= np .abs (
226
+ ss .norm .ppf (0.025 , y_pred , np .sqrt (MSE ))) ).all ()) #check that true value is within 95% conf. int.
227
+
228
+
0 commit comments