Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/descriptivestatistics.cs @ 2806

Last change on this file since 2806 was 2806, checked in by gkronber, 14 years ago

Added plugin for new version of ALGLIB. #875 (Update ALGLIB sources)

File size: 12.4 KB
Line 
1/*************************************************************************
2Copyright (c) 2007, Sergey Bochkanov (ALGLIB project).
3
4>>> SOURCE LICENSE >>>
5This program is free software; you can redistribute it and/or modify
6it under the terms of the GNU General Public License as published by
7the Free Software Foundation (www.fsf.org); either version 2 of the
8License, or (at your option) any later version.
9
10This program is distributed in the hope that it will be useful,
11but WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13GNU General Public License for more details.
14
15A copy of the GNU General Public License is available at
16http://www.fsf.org/licensing/licenses
17
18>>> END OF LICENSE >>>
19*************************************************************************/
20
21using System;
22
23namespace alglib
24{
25    public class descriptivestatistics
26    {
27        /*************************************************************************
28        Calculation of the distribution moments: mean, variance, slewness, kurtosis.
29
30        Input parameters:
31            X       -   sample. Array with whose indexes range within [0..N-1]
32            N       -   sample size.
33           
34        Output parameters:
35            Mean    -   mean.
36            Variance-   variance.
37            Skewness-   skewness (if variance<>0; zero otherwise).
38            Kurtosis-   kurtosis (if variance<>0; zero otherwise).
39
40          -- ALGLIB --
41             Copyright 06.09.2006 by Bochkanov Sergey
42        *************************************************************************/
43        public static void calculatemoments(ref double[] x,
44            int n,
45            ref double mean,
46            ref double variance,
47            ref double skewness,
48            ref double kurtosis)
49        {
50            int i = 0;
51            double v = 0;
52            double v1 = 0;
53            double v2 = 0;
54            double stddev = 0;
55
56            mean = 0;
57            variance = 0;
58            skewness = 0;
59            kurtosis = 0;
60            stddev = 0;
61            if( n<=0 )
62            {
63                return;
64            }
65           
66            //
67            // Mean
68            //
69            for(i=0; i<=n-1; i++)
70            {
71                mean = mean+x[i];
72            }
73            mean = mean/n;
74           
75            //
76            // Variance (using corrected two-pass algorithm)
77            //
78            if( n!=1 )
79            {
80                v1 = 0;
81                for(i=0; i<=n-1; i++)
82                {
83                    v1 = v1+AP.Math.Sqr(x[i]-mean);
84                }
85                v2 = 0;
86                for(i=0; i<=n-1; i++)
87                {
88                    v2 = v2+(x[i]-mean);
89                }
90                v2 = AP.Math.Sqr(v2)/n;
91                variance = (v1-v2)/(n-1);
92                if( (double)(variance)<(double)(0) )
93                {
94                    variance = 0;
95                }
96                stddev = Math.Sqrt(variance);
97            }
98           
99            //
100            // Skewness and kurtosis
101            //
102            if( (double)(stddev)!=(double)(0) )
103            {
104                for(i=0; i<=n-1; i++)
105                {
106                    v = (x[i]-mean)/stddev;
107                    v2 = AP.Math.Sqr(v);
108                    skewness = skewness+v2*v;
109                    kurtosis = kurtosis+AP.Math.Sqr(v2);
110                }
111                skewness = skewness/n;
112                kurtosis = kurtosis/n-3;
113            }
114        }
115
116
117        /*************************************************************************
118        ADev
119
120        Input parameters:
121            X   -   sample (array indexes: [0..N-1])
122            N   -   sample size
123           
124        Output parameters:
125            ADev-   ADev
126
127          -- ALGLIB --
128             Copyright 06.09.2006 by Bochkanov Sergey
129        *************************************************************************/
130        public static void calculateadev(ref double[] x,
131            int n,
132            ref double adev)
133        {
134            int i = 0;
135            double mean = 0;
136
137            mean = 0;
138            adev = 0;
139            if( n<=0 )
140            {
141                return;
142            }
143           
144            //
145            // Mean
146            //
147            for(i=0; i<=n-1; i++)
148            {
149                mean = mean+x[i];
150            }
151            mean = mean/n;
152           
153            //
154            // ADev
155            //
156            for(i=0; i<=n-1; i++)
157            {
158                adev = adev+Math.Abs(x[i]-mean);
159            }
160            adev = adev/n;
161        }
162
163
164        /*************************************************************************
165        Median calculation.
166
167        Input parameters:
168            X   -   sample (array indexes: [0..N-1])
169            N   -   sample size
170
171        Output parameters:
172            Median
173
174          -- ALGLIB --
175             Copyright 06.09.2006 by Bochkanov Sergey
176        *************************************************************************/
177        public static void calculatemedian(double[] x,
178            int n,
179            ref double median)
180        {
181            int i = 0;
182            int ir = 0;
183            int j = 0;
184            int l = 0;
185            int midp = 0;
186            int k = 0;
187            double a = 0;
188            double tval = 0;
189
190            x = (double[])x.Clone();
191
192           
193            //
194            // Some degenerate cases
195            //
196            median = 0;
197            if( n<=0 )
198            {
199                return;
200            }
201            if( n==1 )
202            {
203                median = x[0];
204                return;
205            }
206            if( n==2 )
207            {
208                median = 0.5*(x[0]+x[1]);
209                return;
210            }
211           
212            //
213            // Common case, N>=3.
214            // Choose X[(N-1)/2]
215            //
216            l = 0;
217            ir = n-1;
218            k = (n-1)/2;
219            while( true )
220            {
221                if( ir<=l+1 )
222                {
223                   
224                    //
225                    // 1 or 2 elements in partition
226                    //
227                    if( ir==l+1 & (double)(x[ir])<(double)(x[l]) )
228                    {
229                        tval = x[l];
230                        x[l] = x[ir];
231                        x[ir] = tval;
232                    }
233                    break;
234                }
235                else
236                {
237                    midp = (l+ir)/2;
238                    tval = x[midp];
239                    x[midp] = x[l+1];
240                    x[l+1] = tval;
241                    if( (double)(x[l])>(double)(x[ir]) )
242                    {
243                        tval = x[l];
244                        x[l] = x[ir];
245                        x[ir] = tval;
246                    }
247                    if( (double)(x[l+1])>(double)(x[ir]) )
248                    {
249                        tval = x[l+1];
250                        x[l+1] = x[ir];
251                        x[ir] = tval;
252                    }
253                    if( (double)(x[l])>(double)(x[l+1]) )
254                    {
255                        tval = x[l];
256                        x[l] = x[l+1];
257                        x[l+1] = tval;
258                    }
259                    i = l+1;
260                    j = ir;
261                    a = x[l+1];
262                    while( true )
263                    {
264                        do
265                        {
266                            i = i+1;
267                        }
268                        while( (double)(x[i])<(double)(a) );
269                        do
270                        {
271                            j = j-1;
272                        }
273                        while( (double)(x[j])>(double)(a) );
274                        if( j<i )
275                        {
276                            break;
277                        }
278                        tval = x[i];
279                        x[i] = x[j];
280                        x[j] = tval;
281                    }
282                    x[l+1] = x[j];
283                    x[j] = a;
284                    if( j>=k )
285                    {
286                        ir = j-1;
287                    }
288                    if( j<=k )
289                    {
290                        l = i;
291                    }
292                }
293            }
294           
295            //
296            // If N is odd, return result
297            //
298            if( n%2==1 )
299            {
300                median = x[k];
301                return;
302            }
303            a = x[n-1];
304            for(i=k+1; i<=n-1; i++)
305            {
306                if( (double)(x[i])<(double)(a) )
307                {
308                    a = x[i];
309                }
310            }
311            median = 0.5*(x[k]+a);
312        }
313
314
315        /*************************************************************************
316        Percentile calculation.
317
318        Input parameters:
319            X   -   sample (array indexes: [0..N-1])
320            N   -   sample size, N>1
321            P   -   percentile (0<=P<=1)
322
323        Output parameters:
324            V   -   percentile
325
326          -- ALGLIB --
327             Copyright 01.03.2008 by Bochkanov Sergey
328        *************************************************************************/
329        public static void calculatepercentile(double[] x,
330            int n,
331            double p,
332            ref double v)
333        {
334            int i1 = 0;
335            double t = 0;
336
337            x = (double[])x.Clone();
338
339            System.Diagnostics.Debug.Assert(n>1, "CalculatePercentile: N<=1!");
340            System.Diagnostics.Debug.Assert((double)(p)>=(double)(0) & (double)(p)<=(double)(1), "CalculatePercentile: incorrect P!");
341            internalstatheapsort(ref x, n);
342            if( (double)(p)==(double)(0) )
343            {
344                v = x[0];
345                return;
346            }
347            if( (double)(p)==(double)(1) )
348            {
349                v = x[n-1];
350                return;
351            }
352            t = p*(n-1);
353            i1 = (int)Math.Floor(t);
354            t = t-(int)Math.Floor(t);
355            v = x[i1]*(1-t)+x[i1+1]*t;
356        }
357
358
359        private static void internalstatheapsort(ref double[] arr,
360            int n)
361        {
362            int i = 0;
363            int k = 0;
364            int t = 0;
365            double tmp = 0;
366
367            if( n==1 )
368            {
369                return;
370            }
371            i = 2;
372            do
373            {
374                t = i;
375                while( t!=1 )
376                {
377                    k = t/2;
378                    if( (double)(arr[k-1])>=(double)(arr[t-1]) )
379                    {
380                        t = 1;
381                    }
382                    else
383                    {
384                        tmp = arr[k-1];
385                        arr[k-1] = arr[t-1];
386                        arr[t-1] = tmp;
387                        t = k;
388                    }
389                }
390                i = i+1;
391            }
392            while( i<=n );
393            i = n-1;
394            do
395            {
396                tmp = arr[i];
397                arr[i] = arr[0];
398                arr[0] = tmp;
399                t = 1;
400                while( t!=0 )
401                {
402                    k = 2*t;
403                    if( k>i )
404                    {
405                        t = 0;
406                    }
407                    else
408                    {
409                        if( k<i )
410                        {
411                            if( (double)(arr[k])>(double)(arr[k-1]) )
412                            {
413                                k = k+1;
414                            }
415                        }
416                        if( (double)(arr[t-1])>=(double)(arr[k-1]) )
417                        {
418                            t = 0;
419                        }
420                        else
421                        {
422                            tmp = arr[k-1];
423                            arr[k-1] = arr[t-1];
424                            arr[t-1] = tmp;
425                            t = k;
426                        }
427                    }
428                }
429                i = i-1;
430            }
431            while( i>=1 );
432        }
433    }
434}
Note: See TracBrowser for help on using the repository browser.