SHOGUN  5.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules
HSIC.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (c) The Shogun Machine Learning Toolbox
3  * Written (w) 2012-2013 Heiko Strathmann
4  * Written (w) 2014 Soumyajit De
5  * All rights reserved.
6  *
7  * Redistribution and use in source and binary forms, with or without
8  * modification, are permitted provided that the following conditions are met:
9  *
10  * 1. Redistributions of source code must retain the above copyright notice, this
11  * list of conditions and the following disclaimer.
12  * 2. Redistributions in binary form must reproduce the above copyright notice,
13  * this list of conditions and the following disclaimer in the documentation
14  * and/or other materials provided with the distribution.
15  *
16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
17  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
18  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
19  * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
20  * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
21  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
22  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
23  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26  *
27  * The views and conclusions contained in the software and documentation are those
28  * of the authors and should not be interpreted as representing official policies,
29  * either expressed or implied, of the Shogun Development Team.
30  */
31 
32 #include <shogun/statistics/HSIC.h>
35 #include <shogun/kernel/Kernel.h>
37 
38 using namespace shogun;
39 
41 {
42  init();
43 }
44 
45 CHSIC::CHSIC(CKernel* kernel_p, CKernel* kernel_q, CFeatures* p,
46  CFeatures* q) :
47  CKernelIndependenceTest(kernel_p, kernel_q, p, q)
48 {
49  init();
50 
51  if (p && q && p->get_num_vectors()!=q->get_num_vectors())
52  {
53  SG_ERROR("Only features with equal number of vectors are currently "
54  "possible\n");
55  }
56  else
57  m_num_features=p->get_num_vectors();
58 }
59 
61 {
62 }
63 
64 void CHSIC::init()
65 {
66  SG_ADD(&m_num_features, "num_features",
67  "Number of features from each of the distributions",
69 
70  m_num_features=0;
71 }
72 
74 {
75  SG_DEBUG("entering!\n");
76 
77  REQUIRE(m_kernel_p && m_kernel_q, "No or only one kernel specified!\n");
78 
79  REQUIRE(m_p && m_q, "features needed!\n")
80 
81  /* compute kernel matrices */
84 
85  /* center matrices (MATLAB: Kc=H*K*H) */
86  K.center();
87 
88  /* compute MATLAB: sum(sum(Kc' .* (L))), which is biased HSIC */
89  index_t m=m_num_features;
90  SG_DEBUG("Number of samples %d!\n", m);
91 
92  float64_t result=0;
93  for (index_t i=0; i<m; ++i)
94  {
95  for (index_t j=0; j<m; ++j)
96  result+=K(j, i)*L(i, j);
97  }
98 
99  /* return m times statistic */
100  result/=m;
101 
102  SG_DEBUG("leaving!\n");
103 
104  return result;
105 }
106 
108 {
109  float64_t result=0;
111  {
112  case HSIC_GAMMA:
113  {
114  /* fit gamma and return cdf at statistic */
116  result=CStatistics::gamma_cdf(statistic, params[0], params[1]);
117  break;
118  }
119 
120  default:
121  /* sampling null is handled there */
122  result=CIndependenceTest::compute_p_value(statistic);
123  break;
124  }
125 
126  return result;
127 }
128 
130 {
131  float64_t result=0;
133  {
134  case HSIC_GAMMA:
135  {
136  /* fit gamma and return inverse cdf at statistic */
138 
139  // alpha, beta are shape and rate parameter
140  result=CStatistics::gamma_inverse_cdf(alpha, params[0], params[1]);
141  break;
142  }
143 
144  default:
145  /* sampling null is handled there */
147  break;
148  }
149 
150  return result;
151 }
152 
154 {
155  REQUIRE(m_kernel_p && m_kernel_q, "No or only one kernel specified!\n");
156 
157  REQUIRE(m_p && m_q, "features needed!\n")
158 
159  index_t m=m_num_features;
160 
161  /* compute kernel matrices */
164 
165  /* compute sum and trace of uncentered kernel matrices, needed later */
166  float64_t trace_K=0;
167  float64_t trace_L=0;
168  float64_t sum_K=0;
169  float64_t sum_L=0;
170  for (index_t i=0; i<m; ++i)
171  {
172  trace_K+=K(i,i);
173  trace_L+=L(i,i);
174  for (index_t j=0; j<m; ++j)
175  {
176  sum_K+=K(i,j);
177  sum_L+=L(i,j);
178  }
179  }
180  SG_DEBUG("sum_K: %f, sum_L: %f, trace_K: %f, trace_L: %f\n", sum_K, sum_L,
181  trace_K, trace_L);
182 
183  /* center both matrices: K=H*K*H, L=H*L*H in MATLAB */
184  K.center();
185  L.center();
186 
187  /* compute the trace of MATLAB: (1/6 * Kc.*Lc).^2 Ü */
188  float64_t trace=0;
189  for (index_t i=0; i<m; ++i)
190  trace+=CMath::pow(K(i,i)*L(i,i), 2);
191 
192  trace/=36.0;
193  SG_DEBUG("trace %f\n", trace)
194 
195  /* compute sum of elements of MATLAB: (1/6 * Kc.*Lc).^2 */
196  float64_t sum=0;
197  for (index_t i=0; i<m; ++i)
198  {
199  for (index_t j=0; j<m; ++j)
200  sum+=CMath::pow(K(i,j)*L(i,j), 2);
201  }
202  sum/=36.0;
203  SG_DEBUG("sum %f\n", sum)
204 
205  /* compute MATLAB: 1/m/(m-1)*(sum(sum(varHSIC)) - sum(diag(varHSIC))),
206  * second term is bias correction */
207  float64_t var_hsic=1.0/m/(m-1)*(sum-trace);
208  SG_DEBUG("1.0/m/(m-1)*(sum-trace): %f\n", var_hsic)
209 
210  /* finally, compute variance of hsic under H0
211  * MATLAB: varHSIC = 72*(m-4)*(m-5)/m/(m-1)/(m-2)/(m-3) * varHSIC */
212  var_hsic=72.0*(m-4)*(m-5)/m/(m-1)/(m-2)/(m-3)*var_hsic;
213  SG_DEBUG("var_hsic: %f\n", var_hsic)
214 
215  /* compute mean of matrices with diagonal elements zero on the base of sums
216  * and trace from K and L which were computed above */
217  float64_t mu_x=1.0/m/(m-1)*(sum_K-trace_K);
218  float64_t mu_y=1.0/m/(m-1)*(sum_L-trace_L);
219  SG_DEBUG("mu_x: %f, mu_y: %f\n", mu_x, mu_y)
220 
221  /* compute mean under H0, MATLAB: 1/m * ( 1 +muX*muY - muX - muY ) */
222  float64_t m_hsic=1.0/m*(1+mu_x*mu_y-mu_x-mu_y);
223  SG_DEBUG("m_hsic: %f\n", m_hsic)
224 
225  /* finally, compute parameters of gamma distirbution */
226  float64_t a=CMath::pow(m_hsic, 2)/var_hsic;
227  float64_t b=var_hsic*m/m_hsic;
228  SG_DEBUG("a: %f, b: %f\n", a, b)
229 
230  SGVector<float64_t> result(2);
231  result[0]=a;
232  result[1]=b;
233 
234  SG_DEBUG("leaving!\n")
235  return result;
236 }
237 
239 {
240  SG_DEBUG("entering!\n")
241 
242  /* replace current kernel via precomputed custom kernel and call superclass
243  * method */
244 
245  /* backup references to old kernels */
246  CKernel* kernel_p=m_kernel_p;
247  CKernel* kernel_q=m_kernel_q;
248 
249  /* init kernels before to be sure that everything is fine
250  * kernel function between two samples from different distributions
251  * is never computed - in fact, they may as well have different features */
252  m_kernel_p->init(m_p, m_p);
253  m_kernel_q->init(m_q, m_q);
254 
255  /* precompute kernel matrices */
256  CCustomKernel* precomputed_p=new CCustomKernel(m_kernel_p);
257  CCustomKernel* precomputed_q=new CCustomKernel(m_kernel_q);
258  SG_REF(precomputed_p);
259  SG_REF(precomputed_q);
260 
261  /* temporarily replace own kernels */
262  m_kernel_p=precomputed_p;
263  m_kernel_q=precomputed_q;
264 
265  /* use superclass sample_null which shuffles the entries for one
266  * distribution using index permutation on rows and columns of
267  * kernel matrix from one distribution, while accessing the other
268  * in its original order and then compute statistic */
270 
271  /* restore kernels */
272  m_kernel_p=kernel_p;
273  m_kernel_q=kernel_q;
274 
275  SG_UNREF(precomputed_p);
276  SG_UNREF(precomputed_q);
277 
278  SG_DEBUG("leaving!\n")
279  return null_samples;
280 }
281 
283 {
285  m_num_features=p->get_num_vectors();
286 }
287 
289 {
291  m_num_features=q->get_num_vectors();
292 }
293 
virtual bool init(CFeatures *lhs, CFeatures *rhs)
Definition: Kernel.cpp:98
int32_t index_t
Definition: common.h:62
virtual float64_t compute_p_value(float64_t statistic)
The Custom Kernel allows for custom user provided kernel matrices.
Definition: CustomKernel.h:36
SGMatrix< float64_t > get_kernel_matrix_L()
virtual SGVector< float64_t > sample_null()
virtual int32_t get_num_vectors() const =0
virtual void set_p(CFeatures *p)
Definition: HSIC.cpp:282
#define SG_ERROR(...)
Definition: SGIO.h:129
#define REQUIRE(x,...)
Definition: SGIO.h:206
static float64_t gamma_cdf(float64_t x, float64_t a, float64_t b)
Definition: Statistics.cpp:395
virtual void set_p(CFeatures *p)
#define SG_REF(x)
Definition: SGObject.h:54
virtual float64_t compute_threshold(float64_t alpha)
virtual void set_q(CFeatures *q)
Definition: HSIC.cpp:288
SGMatrix< float64_t > get_kernel_matrix_K()
virtual SGVector< float64_t > sample_null()
Definition: HSIC.cpp:238
SGVector< float64_t > fit_null_gamma()
Definition: HSIC.cpp:153
double float64_t
Definition: common.h:50
virtual float64_t compute_p_value(float64_t statistic)
Definition: HSIC.cpp:107
static float64_t gamma_inverse_cdf(float64_t p, float64_t a, float64_t b)
Definition: Statistics.cpp:513
virtual float64_t compute_threshold(float64_t alpha)
Definition: HSIC.cpp:129
#define SG_UNREF(x)
Definition: SGObject.h:55
#define SG_DEBUG(...)
Definition: SGIO.h:107
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
virtual ~CHSIC()
Definition: HSIC.cpp:60
The class Features is the base class of all feature objects.
Definition: Features.h:68
ENullApproximationMethod m_null_approximation_method
Kernel independence test base class. Provides an interface for performing an independence test...
The Kernel base class.
Definition: Kernel.h:159
#define SG_ADD(...)
Definition: SGObject.h:84
virtual float64_t compute_statistic()
Definition: HSIC.cpp:73
static int32_t pow(bool x, int32_t n)
Definition: Math.h:535
virtual void set_q(CFeatures *q)

SHOGUN Machine Learning Toolbox - Documentation