SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
QuadraticTimeMMD.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 
36 #include <shogun/kernel/Kernel.h>
39 
40 using namespace shogun;
41 
42 #ifdef HAVE_EIGEN3
44 
45 using namespace Eigen;
46 #endif
47 
49 {
50  init();
51 }
52 
54  index_t m) :
55  CKernelTwoSampleTest(kernel, p_and_q, m)
56 {
57  init();
58 }
59 
61  CFeatures* q) : CKernelTwoSampleTest(kernel, p, q)
62 {
63  init();
64 }
65 
67  CKernelTwoSampleTest(custom_kernel, NULL, m)
68 {
69  init();
70 }
71 
73 {
74 }
75 
76 void CQuadraticTimeMMD::init()
77 {
78  SG_ADD(&m_num_samples_spectrum, "num_samples_spectrum", "Number of samples"
79  " for spectrum method null-distribution approximation",
81  SG_ADD(&m_num_eigenvalues_spectrum, "num_eigenvalues_spectrum", "Number of "
82  " Eigenvalues for spectrum method null-distribution approximation",
84  SG_ADD((machine_int_t*)&m_statistic_type, "statistic_type",
85  "Biased or unbiased MMD", MS_NOT_AVAILABLE);
86 
90 }
91 
93  int m, int n)
94 {
95  SG_DEBUG("Entering!\n");
96 
97  /* init kernel with features. NULL check is handled in compute_statistic */
99 
100  /* computing kernel values and their sums on the fly that are used both in
101  computing statistic and variance */
102 
103  /* the following matrix stores row-wise sum of kernel values k(X,X') in
104  the first column and row-wise squared sum of kernel values k^2(X,X')
105  in the second column. m entries in both column */
106  SGMatrix<float64_t> xx_sum_sq_sum_rowwise=m_kernel->
107  row_wise_sum_squared_sum_symmetric_block(0, m);
108 
109  /* row-wise sum of kernel values k(Y,Y'), n entries */
110  SGVector<float64_t> yy_sum_rowwise=m_kernel->
111  row_wise_sum_symmetric_block(m, n);
112 
113  /* row-wise and col-wise sum of kernel values k(X,Y), m+n entries
114  first m entries are row-wise sum, rest n entries are col-wise sum */
115  SGVector<float64_t> xy_sum_rowcolwise=m_kernel->
116  row_col_wise_sum_block(0, m, m, n);
117 
118  /* computing overall sum and squared sum from above for convenience */
119 
120  SGVector<float64_t> xx_sum_rowwise(m);
121  std::copy(xx_sum_sq_sum_rowwise.matrix, xx_sum_sq_sum_rowwise.matrix+m,
122  xx_sum_rowwise.vector);
123 
124  SGVector<float64_t> xy_sum_rowwise(m);
125  std::copy(xy_sum_rowcolwise.vector, xy_sum_rowcolwise.vector+m,
126  xy_sum_rowwise.vector);
127 
128  SGVector<float64_t> xy_sum_colwise(n);
129  std::copy(xy_sum_rowcolwise.vector+m, xy_sum_rowcolwise.vector+m+n,
130  xy_sum_colwise.vector);
131 
132  float64_t xx_sq_sum=0.0;
133  for (index_t i=0; i<xx_sum_sq_sum_rowwise.num_rows; i++)
134  xx_sq_sum+=xx_sum_sq_sum_rowwise(i, 1);
135 
136  float64_t xx_sum=0.0;
137  for (index_t i=0; i<xx_sum_rowwise.vlen; i++)
138  xx_sum+=xx_sum_rowwise[i];
139 
140  float64_t yy_sum=0.0;
141  for (index_t i=0; i<yy_sum_rowwise.vlen; i++)
142  yy_sum+=yy_sum_rowwise[i];
143 
144  float64_t xy_sum=0.0;
145  for (index_t i=0; i<xy_sum_rowwise.vlen; i++)
146  xy_sum+=xy_sum_rowwise[i];
147 
148  /* compute statistic */
149 
150  /* split computations into three terms from JLMR paper (see documentation )*/
151 
152  /* first term */
153  float64_t first=xx_sum/m/(m-1);
154 
155  /* second term */
156  float64_t second=yy_sum/n/(n-1);
157 
158  /* third term */
159  float64_t third=2.0*xy_sum/m/n;
160 
161  float64_t statistic=first+second-third;
162 
163  SG_INFO("Computed statistic!\n");
164  SG_DEBUG("first=%f, second=%f, third=%f\n", first, second, third);
165 
166  /* compute variance under null */
167 
168  /* split computations into three terms (see documentation) */
169 
170  /* first term */
171  float64_t kappa_0=CMath::sq(xx_sum/m/(m-1));
172 
173  /* second term */
174  float64_t kappa_1=0.0;
175  for (index_t i=0; i<m; ++i)
176  kappa_1+=CMath::sq(xx_sum_rowwise[i]/(m-1));
177  kappa_1/=m;
178 
179  /* third term */
180  float64_t kappa_2=xx_sq_sum/m/(m-1);
181 
182  float64_t var_null=2*(kappa_0-2*kappa_1+kappa_2);
183 
184  SG_INFO("Computed variance under null!\n");
185  SG_DEBUG("kappa_0=%f, kappa_1=%f, kappa_2=%f\n", kappa_0, kappa_1, kappa_2);
186 
187  /* compute variance under alternative */
188 
189  /* split computations into four terms (see documentation) */
190 
191  /* first term */
192  float64_t alt_var_first=0.0;
193  for (index_t i=0; i<m; ++i)
194  {
195  // use row-wise sum from k(X,X') and k(X,Y) blocks
196  float64_t term=xx_sum_rowwise[i]/(m-1)-xy_sum_rowwise[i]/n;
197  alt_var_first+=CMath::sq(term);
198  }
199  alt_var_first/=m;
200 
201  /* second term */
202  float64_t alt_var_second=CMath::sq(xx_sum/m/(m-1)-xy_sum/m/n);
203 
204  /* third term */
205  float64_t alt_var_third=0.0;
206  for (index_t i=0; i<n; ++i)
207  {
208  // use row-wise sum from k(Y,Y') and col-wise sum from k(X,Y)
209  // blocks to simulate row-wise sum from k(Y,X) blocks
210  float64_t term=yy_sum_rowwise[i]/(n-1)-xy_sum_colwise[i]/m;
211  alt_var_third+=CMath::sq(term);
212  }
213  alt_var_third/=n;
214 
215  /* fourth term */
216  float64_t alt_var_fourth=CMath::sq(yy_sum/n/(n-1)-xy_sum/m/n);
217 
218  /* finally computing variance */
219  float64_t rho_x=float64_t(m)/(m+n);
220  float64_t rho_y=float64_t(n)/(m+n);
221 
222  float64_t var_alt=4*rho_x*(alt_var_first-alt_var_second)+
223  4*rho_y*(alt_var_third-alt_var_fourth);
224 
225  SG_INFO("Computed variance under alternative!\n");
226  SG_DEBUG("first=%f, second=%f, third=%f, fourth=%f\n", alt_var_first,
227  alt_var_second, alt_var_third, alt_var_fourth);
228 
229  SGVector<float64_t> results(3);
230  results[0]=statistic;
231  results[1]=var_null;
232  results[2]=var_alt;
233 
234  SG_DEBUG("Leaving!\n");
235 
236  return results;
237 }
238 
240 {
241  SG_DEBUG("Entering!\n");
242 
243  /* init kernel with features. NULL check is handled in compute_statistic */
245 
246  /* computing kernel values and their sums on the fly that are used both in
247  computing statistic and variance */
248 
249  /* the following matrix stores row-wise sum of kernel values k(X,X') in
250  the first column and row-wise squared sum of kernel values k^2(X,X')
251  in the second column. m entries in both column */
252  SGMatrix<float64_t> xx_sum_sq_sum_rowwise=m_kernel->
253  row_wise_sum_squared_sum_symmetric_block(0, m, false);
254 
255  /* row-wise sum of kernel values k(Y,Y'), n entries */
256  SGVector<float64_t> yy_sum_rowwise=m_kernel->
257  row_wise_sum_symmetric_block(m, n, false);
258 
259  /* row-wise and col-wise sum of kernel values k(X,Y), m+n entries
260  first m entries are row-wise sum, rest n entries are col-wise sum */
261  SGVector<float64_t> xy_sum_rowcolwise=m_kernel->
262  row_col_wise_sum_block(0, m, m, n);
263 
264  /* computing overall sum and squared sum from above for convenience */
265 
266  SGVector<float64_t> xx_sum_rowwise(m);
267  std::copy(xx_sum_sq_sum_rowwise.matrix, xx_sum_sq_sum_rowwise.matrix+m,
268  xx_sum_rowwise.vector);
269 
270  SGVector<float64_t> xy_sum_rowwise(m);
271  std::copy(xy_sum_rowcolwise.vector, xy_sum_rowcolwise.vector+m,
272  xy_sum_rowwise.vector);
273 
274  SGVector<float64_t> xy_sum_colwise(n);
275  std::copy(xy_sum_rowcolwise.vector+m, xy_sum_rowcolwise.vector+m+n,
276  xy_sum_colwise.vector);
277 
278  float64_t xx_sq_sum=0.0;
279  for (index_t i=0; i<xx_sum_sq_sum_rowwise.num_rows; i++)
280  xx_sq_sum+=xx_sum_sq_sum_rowwise(i, 1);
281 
282  float64_t xx_sum=0.0;
283  for (index_t i=0; i<xx_sum_rowwise.vlen; i++)
284  xx_sum+=xx_sum_rowwise[i];
285 
286  float64_t yy_sum=0.0;
287  for (index_t i=0; i<yy_sum_rowwise.vlen; i++)
288  yy_sum+=yy_sum_rowwise[i];
289 
290  float64_t xy_sum=0.0;
291  for (index_t i=0; i<xy_sum_rowwise.vlen; i++)
292  xy_sum+=xy_sum_rowwise[i];
293 
294  /* compute statistic */
295 
296  /* split computations into three terms from JLMR paper (see documentation )*/
297 
298  /* first term */
299  float64_t first=xx_sum/m/m;
300 
301  /* second term */
302  float64_t second=yy_sum/n/n;
303 
304  /* third term */
305  float64_t third=2.0*xy_sum/m/n;
306 
307  float64_t statistic=first+second-third;
308 
309  SG_INFO("Computed statistic!\n");
310  SG_DEBUG("first=%f, second=%f, third=%f\n", first, second, third);
311 
312  /* compute variance under null */
313 
314  /* split computations into three terms (see documentation) */
315 
316  /* first term */
317  float64_t kappa_0=CMath::sq(xx_sum/m/m);
318 
319  /* second term */
320  float64_t kappa_1=0.0;
321  for (index_t i=0; i<m; ++i)
322  kappa_1+=CMath::sq(xx_sum_rowwise[i]/m);
323  kappa_1/=m;
324 
325  /* third term */
326  float64_t kappa_2=xx_sq_sum/m/m;
327 
328  float64_t var_null=2*(kappa_0-2*kappa_1+kappa_2);
329 
330  SG_INFO("Computed variance under null!\n");
331  SG_DEBUG("kappa_0=%f, kappa_1=%f, kappa_2=%f\n", kappa_0, kappa_1, kappa_2);
332 
333  /* compute variance under alternative */
334 
335  /* split computations into four terms (see documentation) */
336 
337  /* first term */
338  float64_t alt_var_first=0.0;
339  for (index_t i=0; i<m; ++i)
340  {
341  // use row-wise sum from k(X,X') and k(X,Y) blocks
342  float64_t term=xx_sum_rowwise[i]/m-xy_sum_rowwise[i]/n;
343  alt_var_first+=CMath::sq(term);
344  }
345  alt_var_first/=m;
346 
347  /* second term */
348  float64_t alt_var_second=CMath::sq(xx_sum/m/m-xy_sum/m/n);
349 
350  /* third term */
351  float64_t alt_var_third=0.0;
352  for (index_t i=0; i<n; ++i)
353  {
354  // use row-wise sum from k(Y,Y') and col-wise sum from k(X,Y)
355  // blocks to simulate row-wise sum from k(Y,X) blocks
356  float64_t term=yy_sum_rowwise[i]/n-xy_sum_colwise[i]/m;
357  alt_var_third+=CMath::sq(term);
358  }
359  alt_var_third/=n;
360 
361  /* fourth term */
362  float64_t alt_var_fourth=CMath::sq(yy_sum/n/n-xy_sum/m/n);
363 
364  /* finally computing variance */
365  float64_t rho_x=float64_t(m)/(m+n);
366  float64_t rho_y=float64_t(n)/(m+n);
367 
368  float64_t var_alt=4*rho_x*(alt_var_first-alt_var_second)+
369  4*rho_y*(alt_var_third-alt_var_fourth);
370 
371  SG_INFO("Computed variance under alternative!\n");
372  SG_DEBUG("first=%f, second=%f, third=%f, fourth=%f\n", alt_var_first,
373  alt_var_second, alt_var_third, alt_var_fourth);
374 
375  SGVector<float64_t> results(3);
376  results[0]=statistic;
377  results[1]=var_null;
378  results[2]=var_alt;
379 
380  SG_DEBUG("Leaving!\n");
381 
382  return results;
383 }
384 
386 {
387  SG_DEBUG("Entering!\n");
388 
389  /* init kernel with features. NULL check is handled in compute_statistic */
391 
392  /* computing kernel values and their sums on the fly that are used both in
393  computing statistic and variance */
394 
395  /* the following matrix stores row-wise sum of kernel values k(X,X') in
396  the first column and row-wise squared sum of kernel values k^2(X,X')
397  in the second column. n entries in both column */
398  SGMatrix<float64_t> xx_sum_sq_sum_rowwise=m_kernel->
399  row_wise_sum_squared_sum_symmetric_block(0, n);
400 
401  /* row-wise sum of kernel values k(Y,Y'), n entries */
402  SGVector<float64_t> yy_sum_rowwise=m_kernel->
403  row_wise_sum_symmetric_block(n, n);
404 
405  /* row-wise and col-wise sum of kernel values k(X,Y), 2n entries
406  first n entries are row-wise sum, rest n entries are col-wise sum */
407  SGVector<float64_t> xy_sum_rowcolwise=m_kernel->
408  row_col_wise_sum_block(0, n, n, n, true);
409 
410  /* computing overall sum and squared sum from above for convenience */
411 
412  SGVector<float64_t> xx_sum_rowwise(n);
413  std::copy(xx_sum_sq_sum_rowwise.matrix, xx_sum_sq_sum_rowwise.matrix+n,
414  xx_sum_rowwise.vector);
415 
416  SGVector<float64_t> xy_sum_rowwise(n);
417  std::copy(xy_sum_rowcolwise.vector, xy_sum_rowcolwise.vector+n,
418  xy_sum_rowwise.vector);
419 
420  SGVector<float64_t> xy_sum_colwise(n);
421  std::copy(xy_sum_rowcolwise.vector+n, xy_sum_rowcolwise.vector+2*n,
422  xy_sum_colwise.vector);
423 
424  float64_t xx_sq_sum=0.0;
425  for (index_t i=0; i<xx_sum_sq_sum_rowwise.num_rows; i++)
426  xx_sq_sum+=xx_sum_sq_sum_rowwise(i, 1);
427 
428  float64_t xx_sum=0.0;
429  for (index_t i=0; i<xx_sum_rowwise.vlen; i++)
430  xx_sum+=xx_sum_rowwise[i];
431 
432  float64_t yy_sum=0.0;
433  for (index_t i=0; i<yy_sum_rowwise.vlen; i++)
434  yy_sum+=yy_sum_rowwise[i];
435 
436  float64_t xy_sum=0.0;
437  for (index_t i=0; i<xy_sum_rowwise.vlen; i++)
438  xy_sum+=xy_sum_rowwise[i];
439 
440  /* compute statistic */
441 
442  /* split computations into three terms from JLMR paper (see documentation )*/
443 
444  /* first term */
445  float64_t first=xx_sum/n/(n-1);
446 
447  /* second term */
448  float64_t second=yy_sum/n/(n-1);
449 
450  /* third term */
451  float64_t third=2.0*xy_sum/n/(n-1);
452 
453  float64_t statistic=first+second-third;
454 
455  SG_INFO("Computed statistic!\n");
456  SG_DEBUG("first=%f, second=%f, third=%f\n", first, second, third);
457 
458  /* compute variance under null */
459 
460  /* split computations into three terms (see documentation) */
461 
462  /* first term */
463  float64_t kappa_0=CMath::sq(xx_sum/n/(n-1));
464 
465  /* second term */
466  float64_t kappa_1=0.0;
467  for (index_t i=0; i<n; ++i)
468  kappa_1+=CMath::sq(xx_sum_rowwise[i]/(n-1));
469  kappa_1/=n;
470 
471  /* third term */
472  float64_t kappa_2=xx_sq_sum/n/(n-1);
473 
474  float64_t var_null=2*(kappa_0-2*kappa_1+kappa_2);
475 
476  SG_INFO("Computed variance under null!\n");
477  SG_DEBUG("kappa_0=%f, kappa_1=%f, kappa_2=%f\n", kappa_0, kappa_1, kappa_2);
478 
479  /* compute variance under alternative */
480 
481  /* split computations into four terms (see documentation) */
482 
483  /* first term */
484  float64_t alt_var_first=0.0;
485  for (index_t i=0; i<n; ++i)
486  {
487  // use row-wise sum from k(X,X') and k(X,Y) blocks
488  float64_t term=(xx_sum_rowwise[i]-xy_sum_rowwise[i])/(n-1);
489  alt_var_first+=CMath::sq(term);
490  }
491  alt_var_first/=n;
492 
493  /* second term */
494  float64_t alt_var_second=CMath::sq(xx_sum/n/(n-1)-xy_sum/n/(n-1));
495 
496  /* third term */
497  float64_t alt_var_third=0.0;
498  for (index_t i=0; i<n; ++i)
499  {
500  // use row-wise sum from k(Y,Y') and col-wise sum from k(X,Y)
501  // blocks to simulate row-wise sum from k(Y,X) blocks
502  float64_t term=(yy_sum_rowwise[i]-xy_sum_colwise[i])/(n-1);
503  alt_var_third+=CMath::sq(term);
504  }
505  alt_var_third/=n;
506 
507  /* fourth term */
508  float64_t alt_var_fourth=CMath::sq(yy_sum/n/(n-1)-xy_sum/n/(n-1));
509 
510  /* finally computing variance */
511  float64_t rho_x=0.5;
512  float64_t rho_y=0.5;
513 
514  float64_t var_alt=4*rho_x*(alt_var_first-alt_var_second)+
515  4*rho_y*(alt_var_third-alt_var_fourth);
516 
517  SG_INFO("Computed variance under alternative!\n");
518  SG_DEBUG("first=%f, second=%f, third=%f, fourth=%f\n", alt_var_first,
519  alt_var_second, alt_var_third, alt_var_fourth);
520 
521  SGVector<float64_t> results(3);
522  results[0]=statistic;
523  results[1]=var_null;
524  results[2]=var_alt;
525 
526  SG_DEBUG("Leaving!\n");
527 
528  return results;
529 }
530 
532 {
533  return compute_unbiased_statistic_variance(m, n)[0];
534 }
535 
537 {
538  return compute_biased_statistic_variance(m, n)[0];
539 }
540 
542 {
544 }
545 
547 {
548  REQUIRE(m_kernel, "No kernel specified!\n")
549 
550  index_t m=m_m;
551  index_t n=0;
552 
553  /* check if kernel is precomputed (custom kernel) */
555  n=m_kernel->get_num_vec_lhs()-m;
556  else
557  {
558  REQUIRE(m_p_and_q, "The samples are not initialized!\n");
559  n=m_p_and_q->get_num_vectors()-m;
560  }
561 
562  SG_DEBUG("Computing MMD with %d samples from p and %d samples from q!\n",
563  m, n);
564 
565  float64_t result=0;
566  switch (m_statistic_type)
567  {
568  case UNBIASED:
569  result=compute_unbiased_statistic(m, n);
570  result*=m*n/float64_t(m+n);
571  break;
572  case UNBIASED_DEPRECATED:
573  result=compute_unbiased_statistic(m, n);
574  result*=m==n ? m : (m+n);
575  break;
576  case BIASED:
577  result=compute_biased_statistic(m, n);
578  result*=m*n/float64_t(m+n);
579  break;
580  case BIASED_DEPRECATED:
581  result=compute_biased_statistic(m, n);
582  result*=m==n? m : (m+n);
583  break;
584  case INCOMPLETE:
585  REQUIRE(m==n, "Only possible with equal number of samples from both"
586  "distribution!\n")
588  result*=n/2;
589  break;
590  default:
591  SG_ERROR("Unknown statistic type!\n");
592  break;
593  }
594 
595  return result;
596 }
597 
599 {
600  REQUIRE(m_kernel, "No kernel specified!\n")
601 
602  index_t m=m_m;
603  index_t n=0;
604 
605  /* check if kernel is precomputed (custom kernel) */
607  n=m_kernel->get_num_vec_lhs()-m;
608  else
609  {
610  REQUIRE(m_p_and_q, "The samples are not initialized!\n");
611  n=m_p_and_q->get_num_vectors()-m;
612  }
613 
614  SG_DEBUG("Computing MMD with %d samples from p and %d samples from q!\n",
615  m, n);
616 
617  SGVector<float64_t> result(2);
618  switch (m_statistic_type)
619  {
620  case UNBIASED:
621  case UNBIASED_DEPRECATED:
622  {
624  result[0]=res[1];
625  result[1]=res[2];
626  break;
627  }
628  case BIASED:
629  case BIASED_DEPRECATED:
630  {
632  result[0]=res[1];
633  result[1]=res[2];
634  break;
635  }
636  case INCOMPLETE:
637  {
638  REQUIRE(m==n, "Only possible with equal number of samples from both"
639  "distribution!\n")
641  result[0]=res[1];
642  result[1]=res[2];
643  break;
644  }
645  default:
646  SG_ERROR("Unknown statistic type!\n");
647  break;
648  }
649 
650  return result;
651 }
652 
654 {
655  return compute_variance()[0];
656 }
657 
659 {
660  return compute_variance()[1];
661 }
662 
664 {
665  SGVector<float64_t> mmds;
666  if (!multiple_kernels)
667  {
668  /* just one mmd result */
669  mmds=SGVector<float64_t>(1);
670  mmds[0]=compute_statistic();
671  }
672  else
673  {
674  REQUIRE(m_kernel, "No kernel specified!\n")
676  "multiple kernels specified, but underlying kernel is not of type "
677  "K_COMBINED\n");
678 
679  /* cast and allocate memory for results */
681  SG_REF(combined);
682  mmds=SGVector<float64_t>(combined->get_num_subkernels());
683 
684  /* iterate through all kernels and compute statistic */
685  /* TODO this might be done in parallel */
686  for (index_t i=0; i<mmds.vlen; ++i)
687  {
688  CKernel* current=combined->get_kernel(i);
689  /* temporarily replace underlying kernel and compute statistic */
690  m_kernel=current;
691  mmds[i]=compute_statistic();
692 
693  SG_UNREF(current);
694  }
695 
696  /* restore combined kernel */
697  m_kernel=combined;
698  SG_UNREF(combined);
699  }
700 
701  return mmds;
702 }
703 
705 {
706  SGMatrix<float64_t> vars;
707  if (!multiple_kernels)
708  {
709  /* just one mmd result */
710  vars=SGMatrix<float64_t>(1, 2);
712  vars(0, 0)=result[0];
713  vars(0, 1)=result[1];
714  }
715  else
716  {
717  REQUIRE(m_kernel, "No kernel specified!\n")
719  "multiple kernels specified, but underlying kernel is not of type "
720  "K_COMBINED\n");
721 
722  /* cast and allocate memory for results */
724  SG_REF(combined);
725  vars=SGMatrix<float64_t>(combined->get_num_subkernels(), 2);
726 
727  /* iterate through all kernels and compute variance */
728  /* TODO this might be done in parallel */
729  for (index_t i=0; i<vars.num_rows; ++i)
730  {
731  CKernel* current=combined->get_kernel(i);
732  /* temporarily replace underlying kernel and compute variance */
733  m_kernel=current;
735  vars(i, 0)=result[0];
736  vars(i, 1)=result[1];
737 
738  SG_UNREF(current);
739  }
740 
741  /* restore combined kernel */
742  m_kernel=combined;
743  SG_UNREF(combined);
744  }
745 
746  return vars;
747 }
748 
750 {
751  float64_t result=0;
752 
754  {
755  case MMD2_SPECTRUM:
756  {
757 #ifdef HAVE_EIGEN3
758  /* get samples from null-distribution and compute p-value of statistic */
761  null_samples.qsort();
762  index_t pos=null_samples.find_position_to_insert(statistic);
763  result=1.0-((float64_t)pos)/null_samples.vlen;
764 #else // HAVE_EIGEN3
765  SG_ERROR("Only possible if shogun is compiled with EIGEN3 enabled\n");
766 #endif // HAVE_EIGEN3
767  break;
768  }
769 
771  {
772 #ifdef HAVE_EIGEN3
773  /* get samples from null-distribution and compute p-value of statistic */
776  null_samples.qsort();
777  index_t pos=null_samples.find_position_to_insert(statistic);
778  result=1.0-((float64_t)pos)/null_samples.vlen;
779 #else // HAVE_EIGEN3
780  SG_ERROR("Only possible if shogun is compiled with EIGEN3 enabled\n");
781 #endif // HAVE_EIGEN3
782  break;
783  }
784 
785  case MMD2_GAMMA:
786  {
787  /* fit gamma and return cdf at statistic */
789  result=CStatistics::gamma_cdf(statistic, params[0], params[1]);
790  break;
791  }
792 
793  default:
794  result=CKernelTwoSampleTest::compute_p_value(statistic);
795  break;
796  }
797 
798  return result;
799 }
800 
802 {
803  float64_t result=0;
804 
806  {
807  case MMD2_SPECTRUM:
808  {
809 #ifdef HAVE_EIGEN3
810  /* get samples from null-distribution and compute threshold */
813  null_samples.qsort();
814  result=null_samples[index_t(CMath::floor(null_samples.vlen*(1-alpha)))];
815 #else // HAVE_EIGEN3
816  SG_ERROR("Only possible if shogun is compiled with EIGEN3 enabled\n");
817 #endif // HAVE_EIGEN3
818  break;
819  }
820 
822  {
823 #ifdef HAVE_EIGEN3
824  /* get samples from null-distribution and compute threshold */
827  null_samples.qsort();
828  result=null_samples[index_t(CMath::floor(null_samples.vlen*(1-alpha)))];
829 #else // HAVE_EIGEN3
830  SG_ERROR("Only possible if shogun is compiled with EIGEN3 enabled\n");
831 #endif // HAVE_EIGEN3
832  break;
833  }
834 
835  case MMD2_GAMMA:
836  {
837  /* fit gamma and return inverse cdf at alpha */
839  result=CStatistics::inverse_gamma_cdf(alpha, params[0], params[1]);
840  break;
841  }
842 
843  default:
844  /* sampling null is handled here */
846  break;
847  }
848 
849  return result;
850 }
851 
852 
853 #ifdef HAVE_EIGEN3
855  index_t num_eigenvalues)
856 {
857  REQUIRE(m_kernel, "(%d, %d): No kernel set!\n", num_samples,
858  num_eigenvalues);
860  "(%d, %d): No features set and no custom kernel in use!\n",
861  num_samples, num_eigenvalues);
862 
863  index_t m=m_m;
864  index_t n=0;
865 
866  /* check if kernel is precomputed (custom kernel) */
868  n=m_kernel->get_num_vec_lhs()-m;
869  else
870  {
871  REQUIRE(m_p_and_q, "The samples are not initialized!\n");
872  n=m_p_and_q->get_num_vectors()-m;
873  }
874 
875  if (num_samples<=2)
876  {
877  SG_ERROR("Number of samples has to be at least 2, "
878  "better in the hundreds");
879  }
880 
881  if (num_eigenvalues>m+n-1)
882  SG_ERROR("Number of Eigenvalues too large\n");
883 
884  if (num_eigenvalues<1)
885  SG_ERROR("Number of Eigenvalues too small\n");
886 
887  /* imaginary matrix K=[K KL; KL' L] (MATLAB notation)
888  * K is matrix for XX, L is matrix for YY, KL is XY, LK is YX
889  * works since X and Y are concatenated here */
892 
893  /* center matrix K=H*K*H */
894  K.center();
895 
896  /* compute eigenvalues and select num_eigenvalues largest ones */
897  Map<MatrixXd> c_kernel_matrix(K.matrix, K.num_rows, K.num_cols);
898  SelfAdjointEigenSolver<MatrixXd> eigen_solver(c_kernel_matrix);
899  REQUIRE(eigen_solver.info()==Eigen::Success,
900  "Eigendecomposition failed!\n");
901  index_t max_num_eigenvalues=eigen_solver.eigenvalues().rows();
902 
903  /* finally, sample from null distribution */
904  SGVector<float64_t> null_samples(num_samples);
905  for (index_t i=0; i<num_samples; ++i)
906  {
907  null_samples[i]=0;
908  for (index_t j=0; j<num_eigenvalues; ++j)
909  {
911 
912  SG_DEBUG("z_j=%f\n", z_j);
913 
914  float64_t multiple=CMath::sq(z_j);
915 
916  /* take largest EV, scale by 1/(m+n) on the fly and take abs value*/
917  float64_t eigenvalue_estimate=CMath::abs(1.0/(m+n)
918  *eigen_solver.eigenvalues()[max_num_eigenvalues-1-j]);
919 
921  multiple-=1;
922 
923  SG_DEBUG("multiple=%f, eigenvalue=%f\n", multiple,
924  eigenvalue_estimate);
925 
926  null_samples[i]+=eigenvalue_estimate*multiple;
927  }
928  }
929 
930  return null_samples;
931 }
932 
934  index_t num_samples, index_t num_eigenvalues)
935 {
936  REQUIRE(m_kernel, "(%d, %d): No kernel set!\n", num_samples,
937  num_eigenvalues);
939  "(%d, %d): No features set and no custom kernel in use!\n",
940  num_samples, num_eigenvalues);
941 
942  index_t m=m_m;
943  index_t n=0;
944 
945  /* check if kernel is precomputed (custom kernel) */
947  n=m_kernel->get_num_vec_lhs()-m;
948  else
949  {
950  REQUIRE(m_p_and_q, "The samples are not initialized!\n");
951  n=m_p_and_q->get_num_vectors()-m;
952  }
953 
954  if (num_samples<=2)
955  {
956  SG_ERROR("Number of samples has to be at least 2, "
957  "better in the hundreds");
958  }
959 
960  if (num_eigenvalues>m+n-1)
961  SG_ERROR("Number of Eigenvalues too large\n");
962 
963  if (num_eigenvalues<1)
964  SG_ERROR("Number of Eigenvalues too small\n");
965 
966  /* imaginary matrix K=[K KL; KL' L] (MATLAB notation)
967  * K is matrix for XX, L is matrix for YY, KL is XY, LK is YX
968  * works since X and Y are concatenated here */
971 
972  /* center matrix K=H*K*H */
973  K.center();
974 
975  /* compute eigenvalues and select num_eigenvalues largest ones */
976  Map<MatrixXd> c_kernel_matrix(K.matrix, K.num_rows, K.num_cols);
977  SelfAdjointEigenSolver<MatrixXd> eigen_solver(c_kernel_matrix);
978  REQUIRE(eigen_solver.info()==Eigen::Success,
979  "Eigendecomposition failed!\n");
980  index_t max_num_eigenvalues=eigen_solver.eigenvalues().rows();
981 
982  /* precomputing terms with rho_x and rho_y of equation 10 in [1]
983  * (see documentation) */
984  float64_t rho_x=float64_t(m)/(m+n);
985  float64_t rho_y=1-rho_x;
986 
987  /* instead of using two Gaussian rv's ~ N(0,1), we'll use just one rv
988  * ~ N(0, 1/rho_x+1/rho_y) (derived from eq 10 in [1]) */
989  float64_t std_dev=CMath::sqrt(1/rho_x+1/rho_y);
990  float64_t inv_rho_x_y=1/(rho_x*rho_y);
991 
992  SG_DEBUG("Using Gaussian samples ~ N(0,%f)\n", std_dev*std_dev);
993 
994  /* finally, sample from null distribution */
995  SGVector<float64_t> null_samples(num_samples);
996  for (index_t i=0; i<num_samples; ++i)
997  {
998  null_samples[i]=0;
999  for (index_t j=0; j<num_eigenvalues; ++j)
1000  {
1001  /* compute the right hand multiple of eq. 10 in [1] using one RV.
1002  * randn_double() gives a sample from N(0,1), we need samples
1003  * from N(0,1/rho_x+1/rho_y) */
1004  float64_t z_j=std_dev*CMath::randn_double();
1005 
1006  SG_DEBUG("z_j=%f\n", z_j);
1007 
1008  float64_t multiple=CMath::pow(z_j, 2);
1009 
1010  /* take largest EV, scale by 1/(m+n) on the fly and take abs value*/
1011  float64_t eigenvalue_estimate=CMath::abs(1.0/(m+n)
1012  *eigen_solver.eigenvalues()[max_num_eigenvalues-1-j]);
1013 
1015  multiple-=inv_rho_x_y;
1016 
1017  SG_DEBUG("multiple=%f, eigenvalue=%f\n", multiple,
1018  eigenvalue_estimate);
1019 
1020  null_samples[i]+=eigenvalue_estimate*multiple;
1021  }
1022  }
1023 
1024  /* when m=n, return m*MMD^2 instead */
1025  if (m==n)
1026  null_samples.scale(0.5);
1027 
1028  return null_samples;
1029 }
1030 #endif // HAVE_EIGEN3
1031 
1033 {
1034  REQUIRE(m_kernel, "No kernel set!\n");
1036  "No features set and no custom kernel in use!\n");
1037 
1038  index_t n=0;
1039 
1040  /* check if kernel is precomputed (custom kernel) */
1043  else
1044  {
1045  REQUIRE(m_p_and_q, "The samples are not initialized!\n");
1047  }
1048  REQUIRE(m_m==n, "Only possible with equal number of samples "
1049  "from both distribution!\n")
1050 
1051  index_t num_data;
1053  num_data=m_kernel->get_num_vec_rhs();
1054  else
1055  num_data=m_p_and_q->get_num_vectors();
1056 
1057  if (m_m!=num_data/2)
1058  SG_ERROR("Currently, only equal sample sizes are supported\n");
1059 
1060  /* evtl. warn user not to use wrong statistic type */
1062  {
1063  SG_WARNING("Note: provided statistic has "
1064  "to be BIASED. Please ensure that! To get rid of warning,"
1065  "call %s::set_statistic_type(BIASED_DEPRECATED)\n", get_name());
1066  }
1067 
1068  /* imaginary matrix K=[K KL; KL' L] (MATLAB notation)
1069  * K is matrix for XX, L is matrix for YY, KL is XY, LK is YX
1070  * works since X and Y are concatenated here */
1072 
1073  /* compute mean under H0 of MMD, which is
1074  * meanMMD = 2/m * ( 1 - 1/m*sum(diag(KL)) );
1075  * in MATLAB.
1076  * Remove diagonals on the fly */
1077  float64_t mean_mmd=0;
1078  for (index_t i=0; i<m_m; ++i)
1079  {
1080  /* virtual KL matrix is in upper right corner of SHOGUN K matrix
1081  * so this sums the diagonal of the matrix between X and Y*/
1082  mean_mmd+=m_kernel->kernel(i, m_m+i);
1083  }
1084  mean_mmd=2.0/m_m*(1.0-1.0/m_m*mean_mmd);
1085 
1086  /* compute variance under H0 of MMD, which is
1087  * varMMD = 2/m/(m-1) * 1/m/(m-1) * sum(sum( (K + L - KL - KL').^2 ));
1088  * in MATLAB, so sum up all elements */
1089  float64_t var_mmd=0;
1090  for (index_t i=0; i<m_m; ++i)
1091  {
1092  for (index_t j=0; j<m_m; ++j)
1093  {
1094  /* dont add diagonal of all pairs of imaginary kernel matrices */
1095  if (i==j || m_m+i==j || m_m+j==i)
1096  continue;
1097 
1098  float64_t to_add=m_kernel->kernel(i, j);
1099  to_add+=m_kernel->kernel(m_m+i, m_m+j);
1100  to_add-=m_kernel->kernel(i, m_m+j);
1101  to_add-=m_kernel->kernel(m_m+i, j);
1102  var_mmd+=CMath::pow(to_add, 2);
1103  }
1104  }
1105  var_mmd*=2.0/m_m/(m_m-1)*1.0/m_m/(m_m-1);
1106 
1107  /* parameters for gamma distribution */
1108  float64_t a=CMath::pow(mean_mmd, 2)/var_mmd;
1109  float64_t b=var_mmd*m_m / mean_mmd;
1110 
1111  SGVector<float64_t> result(2);
1112  result[0]=a;
1113  result[1]=b;
1114 
1115  return result;
1116 }
1117 
1119  num_samples_spectrum)
1120 {
1121  m_num_samples_spectrum=num_samples_spectrum;
1122 }
1123 
1125  index_t num_eigenvalues_spectrum)
1126 {
1127  m_num_eigenvalues_spectrum=num_eigenvalues_spectrum;
1128 }
1129 
1131  statistic_type)
1132 {
1133  m_statistic_type=statistic_type;
1134 }
1135 

SHOGUN Machine Learning Toolbox - Documentation