- メインページ
- モジュール
- ネームスペース
- データ構造
- ファイル

- ファイル一覧
- グローバル

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