iCub-main
outliersDetection.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2006-2018 Istituto Italiano di Tecnologia (IIT)
3  * Copyright (C) 2006-2010 RobotCub Consortium
4  * All rights reserved.
5  *
6  * This software may be modified and distributed under the terms
7  * of the BSD-3-Clause license. See the accompanying LICENSE file for
8  * details.
9 */
10 
11 #include <cmath>
13 
14 using namespace std;
15 using namespace yarp::os;
16 using namespace yarp::sig;
17 using namespace iCub::ctrl;
18 
19 
20 /**********************************************************************/
21 ModifiedThompsonTau::ModifiedThompsonTau()
22 {
23  // MATLAB code to find out tau values:
24  // n: degrees of freedom
25  // t: the critical student's t value
26  // t=tinv(1-0.05/2,n-2);
27  // tau=(t*(n-1))/sqrt(n*(n-2+t^2));
28 
29  tauLUP[3] =1.1511;
30  tauLUP[4] =1.4250;
31  tauLUP[5] =1.5712;
32  tauLUP[6] =1.6563;
33  tauLUP[7] =1.7110;
34  tauLUP[8] =1.7491;
35  tauLUP[9] =1.7770;
36  tauLUP[10] =1.7984;
37  tauLUP[11] =1.8153;
38  tauLUP[12] =1.8290;
39  tauLUP[13] =1.8403;
40  tauLUP[14] =1.8498;
41  tauLUP[15] =1.8579;
42  tauLUP[16] =1.8649;
43  tauLUP[17] =1.8710;
44  tauLUP[18] =1.8764;
45  tauLUP[19] =1.8811;
46  tauLUP[20] =1.8853;
47  tauLUP[21] =1.8891;
48  tauLUP[22] =1.8926;
49  tauLUP[23] =1.8957;
50  tauLUP[24] =1.8985;
51  tauLUP[25] =1.9011;
52  tauLUP[26] =1.9035;
53  tauLUP[27] =1.9057;
54  tauLUP[28] =1.9078;
55  tauLUP[29] =1.9096;
56  tauLUP[30] =1.9114;
57  tauLUP[31] =1.9130;
58  tauLUP[32] =1.9146;
59  tauLUP[33] =1.9160;
60  tauLUP[34] =1.9174;
61  tauLUP[35] =1.9186;
62  tauLUP[36] =1.9198;
63  tauLUP[37] =1.9209;
64  tauLUP[38] =1.9220;
65  tauLUP[39] =1.9230;
66  tauLUP[40] =1.9240;
67  tauLUP[41] =1.9248;
68  tauLUP[42] =1.9257;
69  tauLUP[43] =1.9265;
70  tauLUP[44] =1.9273;
71  tauLUP[45] =1.9280;
72  tauLUP[46] =1.9288;
73  tauLUP[47] =1.9294;
74  tauLUP[48] =1.9301;
75  tauLUP[49] =1.9307;
76  tauLUP[50] =1.9313;
77  tauLUP[51] =1.9319;
78  tauLUP[52] =1.9324;
79  tauLUP[53] =1.9330;
80  tauLUP[54] =1.9335;
81  tauLUP[55] =1.9340;
82  tauLUP[56] =1.9345;
83  tauLUP[57] =1.9349;
84  tauLUP[58] =1.9354;
85  tauLUP[59] =1.9358;
86  tauLUP[60] =1.9362;
87  tauLUP[61] =1.9366;
88  tauLUP[62] =1.9370;
89  tauLUP[63] =1.9373;
90  tauLUP[64] =1.9377;
91  tauLUP[65] =1.9381;
92  tauLUP[66] =1.9384;
93  tauLUP[67] =1.9387;
94  tauLUP[68] =1.9390;
95  tauLUP[69] =1.9393;
96  tauLUP[70] =1.9396;
97  tauLUP[71] =1.9399;
98  tauLUP[72] =1.9402;
99  tauLUP[73] =1.9405;
100  tauLUP[74] =1.9408;
101  tauLUP[75] =1.9410;
102  tauLUP[76] =1.9413;
103  tauLUP[77] =1.9415;
104  tauLUP[78] =1.9418;
105  tauLUP[79] =1.9420;
106  tauLUP[80] =1.9422;
107  tauLUP[81] =1.9424;
108  tauLUP[82] =1.9427;
109  tauLUP[83] =1.9429;
110  tauLUP[84] =1.9431;
111  tauLUP[85] =1.9433;
112  tauLUP[86] =1.9435;
113  tauLUP[87] =1.9437;
114  tauLUP[88] =1.9439;
115  tauLUP[89] =1.9440;
116  tauLUP[90] =1.9442;
117  tauLUP[91] =1.9444;
118  tauLUP[92] =1.9446;
119  tauLUP[93] =1.9447;
120  tauLUP[94] =1.9449;
121  tauLUP[95] =1.9451;
122  tauLUP[96] =1.9452;
123  tauLUP[97] =1.9454;
124  tauLUP[98] =1.9455;
125  tauLUP[99] =1.9457;
126  tauLUP[100] =1.9458;
127  tauLUP[150] =1.9506;
128  tauLUP[200] =1.9530;
129  tauLUP[500] =1.9572;
130  tauLUP[1000]=1.9586;
131  tauLUP[5000]=1.9597;
132  // n -> inf => tau -> 1.96
133 }
134 
135 
136 /**********************************************************************/
137 set<size_t> ModifiedThompsonTau::detect(const Vector &data, const Property &options)
138 {
139  set<size_t> res;
140  if (data.length()<3)
141  return res;
142 
143  double mean=0.0;
144  double stdev=0.0;
145  size_t N=data.length()-recurIdx.size();
146 
147  if (options.check("mean") && options.check("std"))
148  {
149  mean=options.find("mean").asFloat64();
150  stdev=options.find("std").asFloat64();
151  }
152  else if (N>0)
153  {
154  // find mean and standard deviation
155  for (size_t i=0; i<data.length(); i++)
156  {
157  // don't account for recursive outliers
158  if (recurIdx.find(i)==recurIdx.end())
159  {
160  mean+=data[i];
161  stdev+=data[i]*data[i];
162  }
163  }
164 
165  mean/=N;
166  stdev=sqrt(stdev/N-mean*mean);
167  }
168 
169  size_t i_check;
170  double delta_check=0.0;
171  if (options.check("sorted"))
172  {
173  // don't account for recursive outliers
174  size_t i_1,i_n;
175  for (i_1=0; recurIdx.find(i_1)!=recurIdx.end(); i_1++);
176  for (i_n=data.length()-1; recurIdx.find(i_n)!=recurIdx.end(); i_n--);
177 
178  double delta_1=fabs(data[i_1]-mean);
179  double delta_n=fabs(data[i_n]-mean);
180 
181  if (delta_1>delta_n)
182  {
183  delta_check=delta_1;
184  i_check=i_1;
185  }
186  else
187  {
188  delta_check=delta_n;
189  i_check=i_n;
190  }
191  }
192  else for (size_t i=0; i<data.length(); i++)
193  {
194  // don't account for recursive outliers
195  if (recurIdx.find(i)==recurIdx.end())
196  {
197  double delta=fabs(data[i]-mean);
198  if (delta>delta_check)
199  {
200  delta_check=delta;
201  i_check=i;
202  }
203  }
204  }
205 
206  // choose tau
207  double tau;
208  if (N>5000)
209  tau=1.96;
210  else if (N>1000)
211  tau=tauLUP[5000];
212  else if (N>500)
213  tau=tauLUP[1000];
214  else if (N>200)
215  tau=tauLUP[500];
216  else if (N>150)
217  tau=tauLUP[200];
218  else if (N>100)
219  tau=tauLUP[150];
220  else
221  tau=tauLUP[N];
222 
223  // perform detection
224  if (delta_check>tau*stdev)
225  {
226  // account for current instance
227  res.insert(i_check);
228 
229  // recursive computation
230  if (options.check("recursive"))
231  {
232  // extend the set of current outliers
233  recurIdx.insert(i_check);
234 
235  // recursive call
236  set<size_t> tmpRes=detect(data,options);
237 
238  // if tmpRes is empty then we're done;
239  // hence clear the temporary variables
240  if (tmpRes.empty())
241  recurIdx.clear();
242  // aggregate results
243  else for (set<size_t>::iterator it=tmpRes.begin(); it!=tmpRes.end(); it++)
244  res.insert(*it);
245  }
246 
247  }
248 
249  return res;
250 }
251 
252 
outliersDetection.h
data
@ data
Definition: ST_M1_dataType.h:64
iCub::ctrl
Definition: adaptWinPolyEstimator.h:37