# source:trunk/sources/ALGLIB/descriptivestatistics.cs@2445

Last change on this file since 2445 was 2445, checked in by gkronber, 13 years ago

Fixed #787 (LinearRegressionOperator uses leastsquares function of ALGLIB instead of linearregression function)

File size: 12.2 KB
Line
1/*************************************************************************
2Copyright (c) 2007, Sergey Bochkanov (ALGLIB project).
3
5This program is free software; you can redistribute it and/or modify
7the Free Software Foundation (www.fsf.org); either version 2 of the
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
17
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( variance<0 )
93                {
94                    variance = 0;
95                }
96                stddev = Math.Sqrt(variance);
97            }
98
99            //
100            // Skewness and kurtosis
101            //
102            if( stddev!=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        /*************************************************************************
119
120        Input parameters:
121            X   -   sample (array indexes: [0..N-1])
122            N   -   sample size
123
124        Output parameters:
126
127          -- ALGLIB --
128             Copyright 06.09.2006 by Bochkanov Sergey
129        *************************************************************************/
130        public static void calculateadev(ref double[] x,
131            int n,
133        {
134            int i = 0;
135            double mean = 0;
136
137            mean = 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            //
155            //
156            for(i=0; i<=n-1; i++)
157            {
159            }
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 temp = 0;
189            double tval = 0;
190
191            x = (double[])x.Clone();
192
193
194            //
195            // Some degenerate cases
196            //
197            median = 0;
198            if( n<=0 )
199            {
200                return;
201            }
202            if( n==1 )
203            {
204                median = x[0];
205                return;
206            }
207            if( n==2 )
208            {
209                median = 0.5*(x[0]+x[1]);
210                return;
211            }
212
213            //
214            // Common case, N>=3.
215            // Choose X[(N-1)/2]
216            //
217            l = 0;
218            ir = n-1;
219            k = (n-1)/2;
220            while( true )
221            {
222                if( ir<=l+1 )
223                {
224
225                    //
226                    // 1 or 2 elements in partition
227                    //
228                    if( ir==l+1 & x[ir]<x[l] )
229                    {
230                        tval = x[l];
231                        x[l] = x[ir];
232                        x[ir] = tval;
233                    }
234                    break;
235                }
236                else
237                {
238                    midp = (l+ir)/2;
239                    tval = x[midp];
240                    x[midp] = x[l+1];
241                    x[l+1] = tval;
242                    if( x[l]>x[ir] )
243                    {
244                        tval = x[l];
245                        x[l] = x[ir];
246                        x[ir] = tval;
247                    }
248                    if( x[l+1]>x[ir] )
249                    {
250                        tval = x[l+1];
251                        x[l+1] = x[ir];
252                        x[ir] = tval;
253                    }
254                    if( x[l]>x[l+1] )
255                    {
256                        tval = x[l];
257                        x[l] = x[l+1];
258                        x[l+1] = tval;
259                    }
260                    i = l+1;
261                    j = ir;
262                    a = x[l+1];
263                    while( true )
264                    {
265                        do
266                        {
267                            i = i+1;
268                        }
269                        while( x[i]<a );
270                        do
271                        {
272                            j = j-1;
273                        }
274                        while( x[j]>a );
275                        if( j<i )
276                        {
277                            break;
278                        }
279                        tval = x[i];
280                        x[i] = x[j];
281                        x[j] = tval;
282                    }
283                    x[l+1] = x[j];
284                    x[j] = a;
285                    if( j>=k )
286                    {
287                        ir = j-1;
288                    }
289                    if( j<=k )
290                    {
291                        l = i;
292                    }
293                }
294            }
295
296            //
297            // If N is odd, return result
298            //
299            if( n%2==1 )
300            {
301                median = x[k];
302                return;
303            }
304            a = x[n-1];
305            for(i=k+1; i<=n-1; i++)
306            {
307                if( x[i]<a )
308                {
309                    a = x[i];
310                }
311            }
312            median = 0.5*(x[k]+a);
313        }
314
315
316        /*************************************************************************
317        Percentile calculation.
318
319        Input parameters:
320            X   -   sample (array indexes: [0..N-1])
321            N   -   sample size, N>1
322            P   -   percentile (0<=P<=1)
323
324        Output parameters:
325            V   -   percentile
326
327          -- ALGLIB --
328             Copyright 01.03.2008 by Bochkanov Sergey
329        *************************************************************************/
330        public static void calculatepercentile(double[] x,
331            int n,
332            double p,
333            ref double v)
334        {
335            int i1 = 0;
336            double t = 0;
337
338            x = (double[])x.Clone();
339
340            System.Diagnostics.Debug.Assert(n>1, "CalculatePercentile: N<=1!");
341            System.Diagnostics.Debug.Assert(p>=0 & p<=1, "CalculatePercentile: incorrect P!");
342            internalstatheapsort(ref x, n);
343            if( p==0 )
344            {
345                v = x[0];
346                return;
347            }
348            if( p==1 )
349            {
350                v = x[n-1];
351                return;
352            }
353            t = p*(n-1);
354            i1 = (int)Math.Floor(t);
355            t = t-(int)Math.Floor(t);
356            v = x[i1]*(1-t)+x[i1+1]*t;
357        }
358
359
360        private static void internalstatheapsort(ref double[] arr,
361            int n)
362        {
363            int i = 0;
364            int j = 0;
365            int k = 0;
366            int t = 0;
367            double tmp = 0;
368
369            if( n==1 )
370            {
371                return;
372            }
373            i = 2;
374            do
375            {
376                t = i;
377                while( t!=1 )
378                {
379                    k = t/2;
380                    if( arr[k-1]>=arr[t-1] )
381                    {
382                        t = 1;
383                    }
384                    else
385                    {
386                        tmp = arr[k-1];
387                        arr[k-1] = arr[t-1];
388                        arr[t-1] = tmp;
389                        t = k;
390                    }
391                }
392                i = i+1;
393            }
394            while( i<=n );
395            i = n-1;
396            do
397            {
398                tmp = arr[i];
399                arr[i] = arr[0];
400                arr[0] = tmp;
401                t = 1;
402                while( t!=0 )
403                {
404                    k = 2*t;
405                    if( k>i )
406                    {
407                        t = 0;
408                    }
409                    else
410                    {
411                        if( k<i )
412                        {
413                            if( arr[k]>arr[k-1] )
414                            {
415                                k = k+1;
416                            }
417                        }
418                        if( arr[t-1]>=arr[k-1] )
419                        {
420                            t = 0;
421                        }
422                        else
423                        {
424                            tmp = arr[k-1];
425                            arr[k-1] = arr[t-1];
426                            arr[t-1] = tmp;
427                            t = k;
428                        }
429                    }
430                }
431                i = i-1;
432            }
433            while( i>=1 );
434        }
435    }
436}
Note: See TracBrowser for help on using the repository browser.