Free cookie consent management tool by TermsFeed Policy Generator

source: branches/Persistence Test/ALGLIB/correlation.cs @ 2515

Last change on this file since 2515 was 2430, checked in by gkronber, 15 years ago

Moved ALGLIB code into a separate plugin. #783

File size: 7.9 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 correlation
26    {
27        /*************************************************************************
28        Pearson product-moment correlation coefficient
29
30        Input parameters:
31            X       -   sample 1 (array indexes: [0..N-1])
32            Y       -   sample 2 (array indexes: [0..N-1])
33            N       -   sample size.
34
35        Result:
36            Pearson product-moment correlation coefficient
37
38          -- ALGLIB --
39             Copyright 09.04.2007 by Bochkanov Sergey
40        *************************************************************************/
41        public static double pearsoncorrelation(ref double[] x,
42            ref double[] y,
43            int n)
44        {
45            double result = 0;
46            int i = 0;
47            double xmean = 0;
48            double ymean = 0;
49            double s = 0;
50            double xv = 0;
51            double yv = 0;
52            double t1 = 0;
53            double t2 = 0;
54
55            xv = 0;
56            yv = 0;
57            if( n<=1 )
58            {
59                result = 0;
60                return result;
61            }
62           
63            //
64            // Mean
65            //
66            xmean = 0;
67            ymean = 0;
68            for(i=0; i<=n-1; i++)
69            {
70                xmean = xmean+x[i];
71                ymean = ymean+y[i];
72            }
73            xmean = xmean/n;
74            ymean = ymean/n;
75           
76            //
77            // numerator and denominator
78            //
79            s = 0;
80            xv = 0;
81            yv = 0;
82            for(i=0; i<=n-1; i++)
83            {
84                t1 = x[i]-xmean;
85                t2 = y[i]-ymean;
86                xv = xv+AP.Math.Sqr(t1);
87                yv = yv+AP.Math.Sqr(t2);
88                s = s+t1*t2;
89            }
90            if( xv==0 | yv==0 )
91            {
92                result = 0;
93            }
94            else
95            {
96                result = s/(Math.Sqrt(xv)*Math.Sqrt(yv));
97            }
98            return result;
99        }
100
101
102        /*************************************************************************
103        Spearman's rank correlation coefficient
104
105        Input parameters:
106            X       -   sample 1 (array indexes: [0..N-1])
107            Y       -   sample 2 (array indexes: [0..N-1])
108            N       -   sample size.
109
110        Result:
111            Spearman's rank correlation coefficient
112
113          -- ALGLIB --
114             Copyright 09.04.2007 by Bochkanov Sergey
115        *************************************************************************/
116        public static double spearmanrankcorrelation(double[] x,
117            double[] y,
118            int n)
119        {
120            double result = 0;
121
122            x = (double[])x.Clone();
123            y = (double[])y.Clone();
124
125            rankx(ref x, n);
126            rankx(ref y, n);
127            result = pearsoncorrelation(ref x, ref y, n);
128            return result;
129        }
130
131
132        /*************************************************************************
133        Internal ranking subroutine
134        *************************************************************************/
135        private static void rankx(ref double[] x,
136            int n)
137        {
138            int i = 0;
139            int j = 0;
140            int k = 0;
141            int t = 0;
142            double tmp = 0;
143            int tmpi = 0;
144            double[] r = new double[0];
145            int[] c = new int[0];
146
147           
148            //
149            // Prepare
150            //
151            if( n<=1 )
152            {
153                x[0] = 1;
154                return;
155            }
156            r = new double[n-1+1];
157            c = new int[n-1+1];
158            for(i=0; i<=n-1; i++)
159            {
160                r[i] = x[i];
161                c[i] = i;
162            }
163           
164            //
165            // sort {R, C}
166            //
167            if( n!=1 )
168            {
169                i = 2;
170                do
171                {
172                    t = i;
173                    while( t!=1 )
174                    {
175                        k = t/2;
176                        if( r[k-1]>=r[t-1] )
177                        {
178                            t = 1;
179                        }
180                        else
181                        {
182                            tmp = r[k-1];
183                            r[k-1] = r[t-1];
184                            r[t-1] = tmp;
185                            tmpi = c[k-1];
186                            c[k-1] = c[t-1];
187                            c[t-1] = tmpi;
188                            t = k;
189                        }
190                    }
191                    i = i+1;
192                }
193                while( i<=n );
194                i = n-1;
195                do
196                {
197                    tmp = r[i];
198                    r[i] = r[0];
199                    r[0] = tmp;
200                    tmpi = c[i];
201                    c[i] = c[0];
202                    c[0] = tmpi;
203                    t = 1;
204                    while( t!=0 )
205                    {
206                        k = 2*t;
207                        if( k>i )
208                        {
209                            t = 0;
210                        }
211                        else
212                        {
213                            if( k<i )
214                            {
215                                if( r[k]>r[k-1] )
216                                {
217                                    k = k+1;
218                                }
219                            }
220                            if( r[t-1]>=r[k-1] )
221                            {
222                                t = 0;
223                            }
224                            else
225                            {
226                                tmp = r[k-1];
227                                r[k-1] = r[t-1];
228                                r[t-1] = tmp;
229                                tmpi = c[k-1];
230                                c[k-1] = c[t-1];
231                                c[t-1] = tmpi;
232                                t = k;
233                            }
234                        }
235                    }
236                    i = i-1;
237                }
238                while( i>=1 );
239            }
240           
241            //
242            // compute tied ranks
243            //
244            i = 0;
245            while( i<=n-1 )
246            {
247                j = i+1;
248                while( j<=n-1 )
249                {
250                    if( r[j]!=r[i] )
251                    {
252                        break;
253                    }
254                    j = j+1;
255                }
256                for(k=i; k<=j-1; k++)
257                {
258                    r[k] = 1+((double)(i+j-1))/(double)(2);
259                }
260                i = j;
261            }
262           
263            //
264            // back to x
265            //
266            for(i=0; i<=n-1; i++)
267            {
268                x[c[i]] = r[i];
269            }
270        }
271    }
272}
Note: See TracBrowser for help on using the repository browser.