statistics.h
説明を見る。
1 //
2 // Copyright (c) 2003-2011, MIST Project, Nagoya University
3 // All rights reserved.
4 //
5 // Redistribution and use in source and binary forms, with or without modification,
6 // are permitted provided that the following conditions are met:
7 //
8 // 1. Redistributions of source code must retain the above copyright notice,
9 // this list of conditions and the following disclaimer.
10 //
11 // 2. Redistributions in binary form must reproduce the above copyright notice,
12 // this list of conditions and the following disclaimer in the documentation
13 // and/or other materials provided with the distribution.
14 //
15 // 3. Neither the name of the Nagoya University nor the names of its contributors
16 // may be used to endorse or promote products derived from this software
17 // without specific prior written permission.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
20 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
21 // FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
22 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
24 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
25 // IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
26 // THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 //
28 
33 #ifndef __INCLUDE_MIST_STATISTICS__
34 #define __INCLUDE_MIST_STATISTICS__
35 
36 #include <iostream>
37 #include <stdexcept>
38 #include <cstdio>
39 
40 #ifndef __INCLUDE_MIST_CONF_H__
41 #include "config/mist_conf.h"
42 #endif
43 
44 #ifndef __INCLUDE_MIST_COLOR_H__
45 #include "config/color.h"
46 #endif
47 
48 #ifndef __INCLUDE_MIST_H__
49 #include "mist.h"
50 #endif
51 
52 #ifndef __INCLUDE_MIST_MATRIX__
53 #include "matrix.h"
54 #endif
55 
56 #ifndef __INCLUDE_MIST_NUMERIC__
57 #include "numeric.h"
58 #endif
59 
60 // mist名前空間の始まり
61 _MIST_BEGIN
62 
63 
64 namespace __utility__
65 {
66  template < bool b >
67  struct ____value_type____
68  {
69  typedef double value_type;
70  };
71 
72  template < >
73  struct ____value_type____< true >
74  {
75  typedef rgb< double > value_type;
76  };
77 
78  template < class T >
79  struct __value_type__
80  {
81  typedef typename ____value_type____< is_color< T >::value >::value_type value_type;
82 
83  static value_type t( const value_type &v ){ return( v ); }
84  };
85 
86  template < class T, class Allocator >
87  struct __value_type__< matrix< T, Allocator > >
88  {
89  typedef matrix< T, Allocator > value_type;
90 
91  static value_type t( const value_type &v ){ return( v.t( ) ); }
92  };
93 }
94 
95 
97 namespace statistics
98 {
106 
113  template < class Array >
114  inline typename __utility__::__value_type__< typename Array::value_type >::value_type average( const Array &a )
115  {
116  typedef typename Array::size_type size_type;
117  typedef typename __utility__::__value_type__< typename Array::value_type >::value_type value_type;
118 
119  if( a.empty( ) )
120  {
121  return( value_type( ) );
122  }
123  else
124  {
125  value_type v = a[ 0 ];
126  for( size_type i = 1 ; i < a.size( ) ; i++ )
127  {
128  v += a[ i ];
129  }
130  return( v / static_cast< double >( a.size( ) ) );
131  }
132  }
133 
134 
142  template < class Array >
143  inline typename __utility__::__value_type__< typename Array::value_type >::value_type variance( const Array &a, const typename __utility__::__value_type__< typename Array::value_type >::value_type &ave )
144  {
145  typedef typename Array::size_type size_type;
146  typedef __utility__::__value_type__< typename Array::value_type > __utility_type__;
147  typedef typename __utility_type__::value_type value_type;
148 
149  if( a.empty( ) )
150  {
151  return( value_type( ) );
152  }
153  else
154  {
155  value_type v;
156  {
157  value_type x = a[ 0 ] - ave;
158  v = x * __utility_type__::t( x );
159  }
160 
161  for( size_type i = 1 ; i < a.size( ) ; i++ )
162  {
163  value_type x = a[ i ] - ave;
164  v += x * __utility_type__::t( x );
165  }
166  return( v / static_cast< double >( a.size( ) ) );
167  }
168  }
169 
170 
177  template < class Array >
178  inline typename __utility__::__value_type__< typename Array::value_type >::value_type variance( const Array &a )
179  {
180  return( variance( a, average( a ) ) );
181  }
182 
183 
197  template < class Array1, class Array2 >
198  bool generate_histogram( const Array1 &in, Array2 &out, typename Array1::value_type min, typename Array1::value_type max, typename Array1::value_type bin )
199  {
200  if( is_same_object( in, out ) || in.empty( ) )
201  {
202  return( false );
203  }
204 
205  typedef typename Array1::size_type size_type;
206  typedef typename Array1::difference_type difference_type;
207  typedef typename Array1::value_type value_type;
208 
209  if( bin <= 0 )
210  {
211  return( false );
212  }
213 
214  if( min > max )
215  {
216  value_type tmp = min;
217  min = max;
218  max = tmp;
219  }
220 
221  size_type i;
222  difference_type num = static_cast< size_type >( ( max - min + 1 ) / bin );
223  out.resize( num );
224 
225  for( i = 0 ; i < out.size( ) ; i++ )
226  {
227  out[ i ] = 0;
228  }
229 
230  for( i = 0 ; i < in.size( ) ; i++ )
231  {
232  difference_type index = static_cast< difference_type >( ( in[ i ] - min ) / bin );
233  if( index >= 0 && index < num )
234  {
235  out[ index ] += 1;
236  }
237  }
238 
239  return( true );
240  }
241 
242 
255  template < class Array1, class Array2 >
256  bool generate_histogram( const Array1 &in, Array2 &out, typename Array1::value_type min, typename Array1::value_type max )
257  {
258  return( generate_histogram( in, out, min, max, 1 ) );
259  }
260 
261 
273  template < class Array1, class Array2 >
274  bool generate_histogram( const Array1 &in, Array2 &out, typename Array1::value_type bin )
275  {
276  typedef typename Array1::size_type size_type;
277  typedef typename Array1::difference_type difference_type;
278  typedef typename Array1::value_type value_type;
279 
280  if( in.empty( ) )
281  {
282  return( false );
283  }
284 
285  value_type min = in[ 0 ];
286  value_type max = in[ 0 ];
287  for( size_type i = 1 ; i < in.size( ) ; i++ )
288  {
289  if( min > in[ i ] )
290  {
291  min = in[ i ];
292  }
293  else if( max < in[ i ] )
294  {
295  max = in[ i ];
296  }
297  }
298 
299  return( generate_histogram( in, out, min, max, bin ) );
300  }
301 
312  template < class Array1, class Array2 >
313  bool generate_histogram( const Array1 &in, Array2 &out )
314  {
315  return( generate_histogram( in, out, 1 ) );
316  }
317 
318 
319 
320 
337  template < class Array1, class Array2, class T, class Allocator >
338  bool generate_histogram(
339  const Array1 &in1, const Array2 &in2, array2< T, Allocator > &out,
340  typename Array1::value_type min1, typename Array1::value_type max1,
341  typename Array1::value_type min2, typename Array1::value_type max2,
342  typename Array1::value_type bin
343  )
344  {
345  if( is_same_object( in1, out ) || is_same_object( in2, out ) || in1.empty( ) || in2.empty( ) || in1.size( ) != in2.size( ) )
346  {
347  return( false );
348  }
349 
350  typedef typename Array1::size_type size_type;
351  typedef typename Array1::difference_type difference_type;
352  typedef typename Array1::value_type value_type;
353 
354  if( bin <= 0 )
355  {
356  return( false );
357  }
358 
359  if( min1 > max1 )
360  {
361  value_type tmp = min1;
362  min1 = max1;
363  max1 = tmp;
364  }
365 
366  if( min2 > max2 )
367  {
368  value_type tmp = min2;
369  min2 = max2;
370  max2 = tmp;
371  }
372 
373  difference_type num1 = static_cast< size_type >( ( max1 - min1 + 1 ) / bin );
374  difference_type num2 = static_cast< size_type >( ( max2 - min2 + 1 ) / bin );
375  out.resize( num1, num2 );
376  out.fill( );
377 
378  for( size_type i = 0 ; i < in1.size( ) ; i++ )
379  {
380  difference_type index1 = static_cast< difference_type >( ( in1[ i ] - min1 ) / bin );
381  difference_type index2 = static_cast< difference_type >( ( in2[ i ] - min2 ) / bin );
382  if( index1 >= 0 && index1 < num1 && index2 >= 0 && index2 < num2 )
383  {
384  out( index1, index2 ) += 1;
385  }
386  }
387 
388  return( true );
389  }
390 
391 
407  template < class Array1, class Array2, class T, class Allocator >
408  bool generate_histogram(
409  const Array1 &in1, const Array2 &in2, array2< T, Allocator > &out,
410  typename Array1::value_type min1, typename Array1::value_type max1,
411  typename Array1::value_type min2, typename Array1::value_type max2
412  )
413  {
414  return( generate_histogram( in1, in2, out, min1, max1, min2, max2, 1 ) );
415  }
416 
417 
430  template < class Array1, class Array2, class T, class Allocator >
431  bool generate_histogram( const Array1 &in1, const Array2 &in2, array2< T, Allocator > &out, typename Array1::value_type bin )
432  {
433  typedef typename Array1::size_type size_type;
434  typedef typename Array1::difference_type difference_type;
435  typedef typename Array1::value_type value_type;
436 
437  if( is_same_object( in1, out ) || is_same_object( in2, out ) || in1.empty( ) || in2.empty( ) || in1.size( ) != in2.size( ) )
438  {
439  return( false );
440  }
441 
442  value_type min1 = in1[ 0 ];
443  value_type max1 = in1[ 0 ];
444  for( size_type i = 1 ; i < in1.size( ) ; i++ )
445  {
446  if( min1 > in1[ i ] )
447  {
448  min1 = in1[ i ];
449  }
450  else if( max1 < in1[ i ] )
451  {
452  max1 = in1[ i ];
453  }
454  }
455 
456  value_type min2 = in2[ 0 ];
457  value_type max2 = in2[ 0 ];
458  for( size_type i = 1 ; i < in2.size( ) ; i++ )
459  {
460  if( min2 > in2[ i ] )
461  {
462  min2 = in2[ i ];
463  }
464  else if( max2 < in2[ i ] )
465  {
466  max2 = in2[ i ];
467  }
468  }
469 
470  return( generate_histogram( in1, in2, out, min1, max1, min2, max2, bin ) );
471  }
472 
473 
485  template < class Array1, class Array2, class T, class Allocator >
486  bool generate_histogram( const Array1 &in1, const Array2 &in2, array2< T, Allocator > &out )
487  {
488  return( generate_histogram( in1, in2, out, 1 ) );
489  }
490 
491 
492 
493  namespace detail
494  {
502  template< typename T, typename Allocator >
503  double central_moment( const array2< T, Allocator > &img, size_t x_order, size_t y_order, double x0, double y0 )
504  {
505  double r = 0;
506  for( size_t y = 0 ; y < img.height() ; ++y )
507  {
508  for( size_t x = 0 ; x < img.width() ; ++x )
509  {
510  r += pow( static_cast< double >( x ) - x0, x_order ) * pow( static_cast< double >( y ) - y0, y_order ) * img( x, y );
511  }
512  }
513  return( r );
514  }
515 
524  template< typename T, typename Allocator >
525  double normalized_central_moment( const array2< T, Allocator > &img, size_t x_order, size_t y_order, double u00, double x0, double y0 )
526  {
527  return( central_moment( img, x_order, y_order, x0, y0 ) / u00 );
528  }
529  }
530 
536  template< typename T, typename Allocator >
537  double moment( const array2< T, Allocator > &img, size_t x_order, size_t y_order )
538  {
539  double r = 0;
540 
541  for( size_t y = 0 ; y < img.height() ; ++y )
542  {
543  for( size_t x = 0 ; x < img.width() ; ++x )
544  {
545  r += pow( x, x_order ) * pow( y, y_order ) * img( x, y );
546  }
547  }
548 
549  return( r );
550  }
551 
557  template< typename T, typename Allocator >
558  double central_moment( const array2< T, Allocator > &img, size_t x_order, size_t y_order )
559  {
560  double u00 = moment( img, 0, 0 );
561  double u10 = moment( img, 1, 0 );
562  double u01 = moment( img, 0, 1 );
563 
564  double x0 = u10 / u00;
565  double y0 = u01 / u00;
566 
567  return( detail::central_moment( img, x_order, y_order, x0, y0 ) );
568  }
569 
570 
571 
577  template< typename T, typename Allocator >
578  double normalized_central_moment( const array2< T, Allocator > &img, size_t x_order, size_t y_order )
579  {
580  double u00 = moment( img, 0, 0 );
581  double u10 = moment( img, 1, 0 );
582  double u01 = moment( img, 0, 1 );
583  double x0 = u10 / u00;
584  double y0 = u01 / u00;
585  return( detail::normalized_central_moment( img, x_order, y_order, u00, x0, y0 ) );
586  }
587 
588 
592  template< typename T, typename Allocator >
593  void hu_moments( const array2< T, Allocator > &img, array1< double > &moments )
594  {
595  double u00 = moment( img, 0, 0 );
596  double u10 = moment( img, 1, 0 );
597  double u01 = moment( img, 0, 1 );
598  double x0 = u10 / u00;
599  double y0 = u01 / u00;
600 
601  // 正規化中心モーメントを計算
602  array2< double > ncm( 4, 4 );
603  for( size_t x = 0 ; x <= 3 ; ++x )
604  {
605  for( size_t y = 0 ; y <= 3 ; ++y )
606  {
607  if( x + y > 3 )
608  {
609  continue;
610  }
611  ncm( x, y ) = detail::normalized_central_moment( img, x, y, u00, x0, y0 );
612  }
613  }
614 
615  // Huモーメントを計算
616  moments.resize( 7 );
617  moments( 0 ) = ncm( 2, 0 ) * ncm( 0, 2 );
618  moments( 1 ) = pow( ncm( 2, 0 ) - ncm( 0, 2 ), 2.0 ) + 4 * pow( ncm( 1, 1 ), 2.0 );
619  moments( 2 ) = pow( ncm( 3, 0 ) - 3 * ncm( 1, 2 ), 2.0 ) + pow( 3 * ncm( 2, 1 ) - ncm( 0, 3 ), 2.0 );
620  moments( 3 ) = pow( ncm( 3, 0 ) * ncm( 1, 2 ), 2.0 ) + pow( ncm( 2, 1 ) + ncm( 0, 3 ), 2.0 );
621  moments( 4 ) = ( ncm( 3, 0 ) - 3 * ncm( 1, 2 ) ) *
622  ( ncm( 3, 0 ) + ncm( 1, 2 ) ) *
623  ( pow( ncm( 3, 0 ) + ncm( 1, 2 ), 2.0 ) - 3 * pow( ncm( 2, 1 ) + ncm( 0, 3 ), 2.0 ) ) +
624  ( 3 * ncm( 2, 1 ) - ncm( 0, 3 ) ) *
625  ( ncm( 2, 1 ) + ncm( 0, 3 ) ) *
626  ( 3 * pow( ncm( 3, 0 ) + ncm( 1, 2 ), 2.0 ) - pow( ncm( 2, 1 ) + ncm( 0, 3 ), 2.0 ) );
627  moments( 5 ) = ( ncm( 2, 0 ) - ncm( 0, 2 ) ) *
628  ( pow( ncm( 3, 0 ) + ncm( 1, 2 ), 2.0 ) - pow( ncm( 2, 1 ) + ncm( 0, 3 ), 2.0 ) ) +
629  4 * ncm( 1, 1 ) * ( ncm( 3, 0 ) + ncm( 1, 2 ) ) * ( ncm( 2, 1 ) + ncm( 0, 3 ) );
630  moments( 6 ) = ( 3 * ncm( 2, 1 ) - ncm( 0, 3 ) ) *
631  ( ncm( 2, 1 ) + ncm( 0, 3 ) ) *
632  ( 3 * pow( ncm( 3, 0 ) + ncm( 1, 2 ), 2.0 ) - pow( ncm( 2, 1 ) + ncm( 0, 3 ), 2.0 ) ) -
633  ( ncm( 3, 0 ) - 3 * ncm( 2, 1 ) ) *
634  ( ncm( 2, 1 ) + ncm( 0, 3 ) ) *
635  ( 3 * pow( ncm( 3, 0 ) + ncm( 1, 2 ), 2.0 ) - pow( ncm( 2, 1 ) + ncm( 0, 3 ), 2.0 ) );
636  }
637 
638  namespace normal_distribution
639  {
645  template< typename T, typename Allocator >
646  double probability( const mist::matrix< T, Allocator > &x, const mist::matrix< T, Allocator > &u, const mist::matrix< T, Allocator > &sigma )
647  {
648  if( x.rows() != u.rows() || u.cols() != 1 || x.rows() != sigma.rows() || sigma.rows() != sigma.cols() )
649  {
650  throw std::invalid_argument( "" );
651  }
652  mist::matrix< T, Allocator > b = x - u;
653  double t = ( b.t() * mist::inverse( sigma ) * b )( 0, 0 );
654  double d = mist::det( sigma );
655  return ( 1.0 / pow( 2 * M_PI, 0.5 * x.rows() ) ) * ( 1.0 / d ) * exp( -0.5 * t );
656  }
657 
662  template< typename T, typename Allocator >
663  void estimate( const mist::matrix< T, Allocator > &sample, mist::matrix< T, Allocator > &average, mist::matrix< T, Allocator > &variance )
664  {
665  if( sample.cols() <= 1 )
666  {
667  throw std::invalid_argument( "" );
668  }
669 
670  average.resize( sample.rows(), 1 );
671  for( size_t i = 0 ; i < sample.cols() ; ++i )
672  {
673  for( size_t j = 0 ; j < sample.rows() ; ++j )
674  {
675  average( j, 0 ) += sample( j, i );
676  }
677  }
678  average /= sample.cols();
679 
680  mist::matrix< T, Allocator > tmp = sample;
681  for( size_t i = 0 ; i < sample.cols() ; ++i )
682  {
683  for( size_t j = 0 ; j < sample.rows() ; ++j )
684  {
685  tmp( j, i ) -= average( j, 0 );
686  }
687  }
688  variance = ( tmp * tmp.t() ) / ( sample.cols() - 1 );
689  }
690  }
691 
692 
693  namespace parzen
694  {
700  template< typename T, typename Allocator >
701  double density( const mist::matrix< T, Allocator > &sample, const mist::matrix< T, Allocator > &x, double band )
702  {
703  if( sample.cols() == 0 )
704  {
705  throw std::invalid_argument( "" );
706  }
707 
708  if( sample.rows() != x.rows() )
709  {
710  throw std::invalid_argument( "" );
711  }
712 
713  if( band <= 0 )
714  {
715  throw std::invalid_argument( "" );
716  }
717 
718  double t = 0.0;
719  for( size_t i = 0 ; i < sample.cols() ; ++i )
720  {
721  double a = 1.0;
722  for( size_t j = 0 ; j < sample.rows() ; ++j )
723  {
724  double s = fabs( sample( j, i ) - x( j, 0 ) ) / band;
725  if( s > 0.5 )
726  {
727  a = 0.0;
728  break;
729  }
730  }
731  t += a;
732  }
733  return ( 1.0 / ( pow( band, sample.rows() ) * sample.cols() ) ) * t;
734  }
735 
741  template< typename T, typename Allocator >
742  double test( const mist::matrix< T, Allocator > &sample, double band, int n = 5 )
743  {
744  if( sample.cols() == 0 )
745  {
746  throw std::invalid_argument( "" );
747  }
748 
749  if( n == 0 )
750  {
751  throw std::invalid_argument( "" );
752  }
753 
754  if( band <= 0 )
755  {
756  throw std::invalid_argument( "" );
757  }
758 
759  int ns = sample.cols() / n;
760  int n1 = ns * ( n - 1 );
761 
762  double J = 0;
763 
764  mist::matrix< T, Allocator > a( sample.rows(), n1 );
765  for( int i = 0 ; i < n ; ++i )
766  {
767  int idx = 0;
768  for( size_t j = 0 ; j < sample.cols() ; ++j )
769  {
770  if( static_cast< int >( j / ns ) != i )
771  {
772  for( size_t k = 0 ; k < sample.rows() ; ++k )
773  {
774  a( k, idx ) = sample( k, j );
775  }
776  ++idx;
777  }
778  }
779 
780  for( size_t j = 0 ; j < sample.cols() ; ++j )
781  {
782  if( static_cast< int >( j / ns ) == i )
783  {
784  mist::matrix< T, Allocator > x( sample.rows(), 1 );
785  sample.trim( x, 0, j, -1, 1 );
786  J += log( density( a, x, band ) );
787  }
788  }
789  }
790  return J / n;
791  }
792  }
793 
794  namespace kernel
795  {
801  template< typename T, typename Allocator >
802  double density( const mist::matrix< T, Allocator > &sample, const mist::matrix< T, Allocator > &x, double b )
803  {
804  double t = 0.0;
805  for( size_t i = 0 ; i < sample.cols() ; ++i )
806  {
807  mist::matrix< T, Allocator > d = x;
808  for( size_t j = 0 ; j < sample.rows() ; ++j )
809  {
810  d( j, 0 ) -= sample( j, i );
811  }
812  d /= b;
813  t += 1.0 / pow( 2.0 * M_PI, 0.5 * sample.rows() ) * exp( -0.5 * ( d.t() * d )( 0, 0 ) );
814  }
815  return ( 1.0 / ( pow( b, sample.rows() ) * sample.cols() ) ) * t;
816  }
817 
823  template< typename T, typename Allocator >
824  double test( const mist::matrix< T, Allocator > &sample, double band, int n = 5 )
825  {
826  if( sample.cols() == 0 )
827  {
828  throw std::invalid_argument( "" );
829  }
830 
831  if( n == 0 )
832  {
833  throw std::invalid_argument( "" );
834  }
835 
836  if( band <= 0 )
837  {
838  throw std::invalid_argument( "" );
839  }
840 
841  int ns = sample.cols() / n;
842  int n1 = ns * ( n - 1 );
843 
844  double J = 0;
845 
846  mist::matrix< T, Allocator > a( sample.rows(), n1 );
847  for( int i = 0 ; i < n ; ++i )
848  {
849  int idx = 0;
850  for( size_t j = 0 ; j < sample.cols() ; ++j )
851  {
852  if( static_cast< int >( j / ns ) != i )
853  {
854  for( size_t k = 0 ; k < sample.rows() ; ++k )
855  {
856  a( k, idx ) = sample( k, j );
857  }
858  ++idx;
859  }
860  }
861 
862  for( size_t j = 0 ; j < sample.cols() ; ++j )
863  {
864  if( static_cast< int >( j / ns ) == i )
865  {
866  mist::matrix< T, Allocator > x( sample.rows(), 1 );
867  sample.trim( x, 0, j, -1, 1 );
868  J += log( density( a, x, band ) );
869  }
870  }
871  }
872  return J / n;
873  }
874  }
875 
876  namespace knn
877  {
883  template< typename T, typename Allocator >
884  double density( const mist::matrix< T, Allocator > &sample, const mist::matrix< T, Allocator > &x, int k )
885  {
886  std::vector< T, Allocator > ssample( sample.size() );
887  for( int i = 0 ; i < static_cast< int >( sample.cols() ) ; ++i )
888  {
889  double dist = 0;
890  for( int j = 0 ; j < static_cast< int >( sample.rows() ) ; ++j )
891  {
892  dist += pow( x( j, 0 ) - sample( j, i ), 2.0 );
893  }
894  ssample[ i ] = dist;
895  }
896 
897  std::sort( ssample.begin(), ssample.end() );
898 
899  double r = sqrt( ssample[ k - 1 ] );
900 
901  int h = 0;
902  int d = static_cast< int >( sample.rows() );
903  double g = 1.0;
904  if( d % 2 == 0 )
905  {
906  d /= 2;
907  while( d != 1 )
908  {
909  g *= d;
910  d -= 1;
911  }
912  }
913  else
914  {
915  g = 0.5;
916  while( d != 1 )
917  {
918  g *= ( 0.5 * d );
919  d -= 2;
920  }
921  h = 1;
922  }
923  return ( g * k ) / ( pow( M_PI, 0.5 * ( sample.rows() - h ) ) * pow( r, sample.rows() ) * sample.cols() );
924  }
925 
931  template< typename T, typename Allocator >
932  double test( const mist::matrix< T, Allocator > &sample, int k, int n = 5 )
933  {
934  if( sample.cols() == 0 )
935  {
936  throw std::invalid_argument( "" );
937  }
938 
939  if( n == 0 )
940  {
941  throw std::invalid_argument( "" );
942  }
943 
944  if( k <= 0 )
945  {
946  throw std::invalid_argument( "" );
947  }
948 
949  int ns = sample.cols() / n;
950  int n1 = ns * ( n - 1 );
951 
952  double J = 0;
953 
954  mist::matrix< T, Allocator > a( sample.rows(), n1 );
955  for( int i = 0 ; i < n ; ++i )
956  {
957  int idx = 0;
958  for( size_t j = 0 ; j < sample.cols() ; ++j )
959  {
960  if( static_cast< int >( j / ns ) != i )
961  {
962  for( size_t k = 0 ; k < sample.rows() ; ++k )
963  {
964  a( k, idx ) = sample( k, j );
965  }
966  ++idx;
967  }
968  }
969 
970  for( size_t j = 0 ; j < sample.cols() ; ++j )
971  {
972  if( static_cast< int >( j / ns ) == i )
973  {
974  mist::matrix< T, Allocator > x( sample.rows(), 1 );
975  sample.trim( x, 0, j, -1, 1 );
976  J += log( density( a, x, k ) );
977  }
978  }
979  }
980  return J / n;
981  }
982  }
984  // 統計処理の終わり
985 }
986 
987 
988 // mist名前空間の終わり
989 _MIST_END
990 
991 
992 #endif // __INCLUDE_MIST_STATISTICS__
993 
Generated on Wed Nov 12 2014 19:44:25 for MIST by doxygen 1.8.1.2