A.I.M.S algorithms


wlms-estimator.h
Go to the documentation of this file.
1 /* This software and supporting documentation are distributed by
2  * Institut Federatif de Recherche 49
3  * CEA/NeuroSpin, Batiment 145,
4  * 91191 Gif-sur-Yvette cedex
5  * France
6  *
7  * This software is governed by the CeCILL-B license under
8  * French law and abiding by the rules of distribution of free software.
9  * You can use, modify and/or redistribute the software under the
10  * terms of the CeCILL-B license as circulated by CEA, CNRS
11  * and INRIA at the following URL "http://www.cecill.info".
12  *
13  * As a counterpart to the access to the source code and rights to copy,
14  * modify and redistribute granted by the license, users are provided only
15  * with a limited warranty and the software's author, the holder of the
16  * economic rights, and the successive licensors have only limited
17  * liability.
18  *
19  * In this respect, the user's attention is drawn to the risks associated
20  * with loading, using, modifying and/or developing or reproducing the
21  * software by the user in light of its specific status of free software,
22  * that may mean that it is complicated to manipulate, and that also
23  * therefore means that it is reserved for developers and experienced
24  * professionals having in-depth computer knowledge. Users are therefore
25  * encouraged to load and test the software's suitability as regards their
26  * requirements in conditions enabling the security of their systems and/or
27  * data to be ensured and, more generally, to use and operate it in the
28  * same conditions as regards security.
29  *
30  * The fact that you are presently reading this means that you have had
31  * knowledge of the CeCILL-B license and that you accept its terms.
32  */
33 
34 
35 #ifndef AIMS_ESTIMATION_WLMS_ESTIMATOR_H
36 #define AIMS_ESTIMATION_WLMS_ESTIMATOR_H
37 
38 #include <cstdlib>
41 #include <aims/estimation/m-func.h>
42 
43 //
44 // Weighted Linear Mean Squared M-estimator
45 //
46 template < int D >
48 {
49  public:
50  WLMSEstimator( MEstimatorFunc& mfunc, float error,
51  bool verbose = false )
52  : MEstimator<D>(), _mfunc( mfunc ), _error( error ),
53  _verbose( verbose) { }
54  virtual ~WLMSEstimator() { }
55 
56  void doit( const AimsData< AimsVector< float, D > >& x,
57  const AimsData< float >& y, float& a,
59 
60  protected:
62  float _error;
63  bool _verbose;
64 
65  float sigma( const AimsData< AimsVector< float, D > >& x,
66  const AimsData< float >& y, float a,
67  const AimsVector< float, D >& b,
68  AimsData<float>& error );
69 
70  void step( const AimsData< AimsVector< float, D > >& x,
71  const AimsData< float >& y, float& a,
73  const AimsData<float>& error );
74 };
75 
76 
77 template < int D > inline
79  const AimsData< float >& y, float a,
80  const AimsVector< float, D >& b,
81  AimsData<float>& error )
82 {
83  int N = y.dimX(), n, k;
84  for ( n = 0; n < N; n++ )
85  {
86  error( n ) = y ( n ) - a;
87  for ( k = 0; k < D; k++ )
88  error( n ) -= b.item( k ) * x( n ).item( k );
89  error( n ) = fabs( error( n ) );
90  }
91 
92  AimsData<float> errorbis = error.clone();
93  incSorting( errorbis );
94 
95  return 1.4826 * errorbis( N / 2 );
96 }
97 
98 
99 template < int D >
100 inline
101 void
103  const AimsData< float >& y, float& a,
105  const AimsData<float>& error )
106 {
107  int N = y.dimX();
108  AimsData<float> mat( D + 1, D + 1 );
109  AimsData<float> vec( D + 1 );
110  AimsData<float> weight( N );
111 
112  mat = 0.0;
113  vec = 0.0;
114  weight = 0.0;
115 
116  int k, n;
117 
118  for ( n = 0; n < N; n++ )
119  {
120  weight( n ) = _mfunc.weight( error( n ) );
121  mat( 0, 0 ) += weight( n );
122  }
123 
124  for ( k = 1; k <= D; k++ )
125  {
126  for ( n = 0; n < N; n++ )
127  mat( k, 0 ) += weight( n ) * x( n ).item( k - 1 );
128  mat( 0, k ) = mat( k, 0 );
129  }
130 
131  int k1, k2;
132  for ( k1 = 1; k1 <= D; k1++ )
133  for ( k2 = 1; k2 <= k1; k2++ )
134  {
135  for ( n = 0; n < N; n++ )
136  mat( k1, k2 ) += weight( n ) *
137  x( n ).item( k1 - 1 ) * x( n ).item( k2 - 1 );
138  mat( k2, k1 ) = mat( k1, k2 );
139  }
140 
141 
142  for ( n = 0; n < N; n++ )
143  {
144  vec( 0 ) += weight( n ) * y( n );
145  for ( k = 1; k <= D; k++ )
146  vec( k ) += weight( n ) * y( n ) * x( n ).item( k - 1 );
147  }
148 
149  AimsData<float> res = AimsLinearResolutionLU( mat, vec );
150  a = res( 0 );
151  for ( k = 1; k <= D; k++ )
152  b.item( k - 1 ) = res( k );
153 }
154 
155 
156 template < int D > inline
158  const AimsData< float >& y, float& a,
160 {
161  ASSERT( x.dimY() == 1 && x.dimZ() == 1 && x.dimT() == 1 );
162  ASSERT( y.dimY() == 1 && y.dimZ() == 1 && y.dimT() == 1 );
163  ASSERT( x.dimX() == y.dimX() );
164 
165  float obja, objb;
166  float old_a;
168 
169  // estimate starting value with LMS fit
170  LMSEstimator<D> lms;
171  lms.doit( x, y, a, b );
172 
173  AimsData<float> error( x.dimX() );
174  int count = 0;
175  // perform iteration of robust estimation until convergence
176  do
177  {
178  ASSERT( count != 100 ); // maximum iteration
179  old_a = a;
180  old_b = b;
181  _mfunc.setSigma( sigma( x, y, old_a,old_b, error ) );
182  step( x, y, a, b, error );
183  obja = fabs( a - old_a ) / ( fabs( a ) + fabs( old_a ) );
184  objb = norm( b - old_b ) / ( norm( b ) + norm( old_b ) );
185  if ( _verbose )
186  std::cout
187  << "Obja : " << obja << " "
188  << "Objb : " << objb << " "
189  << "a : " << a << " "
190  << "b : " << b << std::endl;
191  count++;
192  }
193  while ( obja > _error || objb > _error );
194 }
195 
196 
197 #endif
virtual void doit(const AimsData< AimsVector< float, D > > &, const AimsData< float > &, float &, AimsVector< float, D > &)
Definition: m-estimator.h:54
AIMSDATA_API AimsData< float > AimsLinearResolutionLU(const AimsData< float > &matrix, const AimsData< float > &b)
MEstimatorFunc & _mfunc
int dimX() const
WLMSEstimator(MEstimatorFunc &mfunc, float error, bool verbose=false)
void doit(const AimsData< AimsVector< float, D > > &x, const AimsData< float > &y, float &a, AimsVector< float, D > &b)
Definition: lms-estimator.h:64
int dimT() const
int verbose
void step(const AimsData< AimsVector< float, D > > &x, const AimsData< float > &y, float &a, AimsVector< float, D > &b, const AimsData< float > &error)
void incSorting(AimsData< T > &thing)
void doit(const AimsData< AimsVector< float, D > > &x, const AimsData< float > &y, float &a, AimsVector< float, D > &b)
AimsData< T > clone() const
const T & item(int d) const
virtual ~WLMSEstimator()
float sigma(const AimsData< AimsVector< float, D > > &x, const AimsData< float > &y, float a, const AimsVector< float, D > &b, AimsData< float > &error)
int dimZ() const
int dimY() const
AIMSDATA_API float norm(const Tensor &thing)
#define AIMSALGO_API
#define ASSERT(EX)