Index: /trunk/sources/ALGLIB/ALGLIB.csproj
===================================================================
--- /trunk/sources/ALGLIB/ALGLIB.csproj	(revision 2444)
+++ /trunk/sources/ALGLIB/ALGLIB.csproj	(revision 2445)
@@ -88,6 +88,11 @@
       <CopyToOutputDirectory>Always</CopyToOutputDirectory>
     </Compile>
+    <Compile Include="descriptivestatistics.cs" />
+    <Compile Include="gammaf.cs" />
+    <Compile Include="igammaf.cs" />
     <Compile Include="leastsquares.cs" />
+    <Compile Include="linreg.cs" />
     <Compile Include="lq.cs" />
+    <Compile Include="normaldistr.cs" />
     <Compile Include="Properties\AssemblyInfo.cs" />
     <Compile Include="qr.cs" />
Index: /trunk/sources/ALGLIB/descriptivestatistics.cs
===================================================================
--- /trunk/sources/ALGLIB/descriptivestatistics.cs	(revision 2445)
+++ /trunk/sources/ALGLIB/descriptivestatistics.cs	(revision 2445)
@@ -0,0 +1,436 @@
+/*************************************************************************
+Copyright (c) 2007, Sergey Bochkanov (ALGLIB project).
+
+>>> SOURCE LICENSE >>>
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation (www.fsf.org); either version 2 of the 
+License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+A copy of the GNU General Public License is available at
+http://www.fsf.org/licensing/licenses
+
+>>> END OF LICENSE >>>
+*************************************************************************/
+
+using System;
+
+namespace alglib
+{
+    public class descriptivestatistics
+    {
+        /*************************************************************************
+        Calculation of the distribution moments: mean, variance, slewness, kurtosis.
+
+        Input parameters:
+            X       -   sample. Array with whose indexes range within [0..N-1]
+            N       -   sample size.
+            
+        Output parameters:
+            Mean    -   mean.
+            Variance-   variance.
+            Skewness-   skewness (if variance<>0; zero otherwise).
+            Kurtosis-   kurtosis (if variance<>0; zero otherwise).
+
+          -- ALGLIB --
+             Copyright 06.09.2006 by Bochkanov Sergey
+        *************************************************************************/
+        public static void calculatemoments(ref double[] x,
+            int n,
+            ref double mean,
+            ref double variance,
+            ref double skewness,
+            ref double kurtosis)
+        {
+            int i = 0;
+            double v = 0;
+            double v1 = 0;
+            double v2 = 0;
+            double stddev = 0;
+
+            mean = 0;
+            variance = 0;
+            skewness = 0;
+            kurtosis = 0;
+            stddev = 0;
+            if( n<=0 )
+            {
+                return;
+            }
+            
+            //
+            // Mean
+            //
+            for(i=0; i<=n-1; i++)
+            {
+                mean = mean+x[i];
+            }
+            mean = mean/n;
+            
+            //
+            // Variance (using corrected two-pass algorithm)
+            //
+            if( n!=1 )
+            {
+                v1 = 0;
+                for(i=0; i<=n-1; i++)
+                {
+                    v1 = v1+AP.Math.Sqr(x[i]-mean);
+                }
+                v2 = 0;
+                for(i=0; i<=n-1; i++)
+                {
+                    v2 = v2+(x[i]-mean);
+                }
+                v2 = AP.Math.Sqr(v2)/n;
+                variance = (v1-v2)/(n-1);
+                if( variance<0 )
+                {
+                    variance = 0;
+                }
+                stddev = Math.Sqrt(variance);
+            }
+            
+            //
+            // Skewness and kurtosis
+            //
+            if( stddev!=0 )
+            {
+                for(i=0; i<=n-1; i++)
+                {
+                    v = (x[i]-mean)/stddev;
+                    v2 = AP.Math.Sqr(v);
+                    skewness = skewness+v2*v;
+                    kurtosis = kurtosis+AP.Math.Sqr(v2);
+                }
+                skewness = skewness/n;
+                kurtosis = kurtosis/n-3;
+            }
+        }
+
+
+        /*************************************************************************
+        ADev
+
+        Input parameters:
+            X   -   sample (array indexes: [0..N-1])
+            N   -   sample size
+            
+        Output parameters:
+            ADev-   ADev
+
+          -- ALGLIB --
+             Copyright 06.09.2006 by Bochkanov Sergey
+        *************************************************************************/
+        public static void calculateadev(ref double[] x,
+            int n,
+            ref double adev)
+        {
+            int i = 0;
+            double mean = 0;
+
+            mean = 0;
+            adev = 0;
+            if( n<=0 )
+            {
+                return;
+            }
+            
+            //
+            // Mean
+            //
+            for(i=0; i<=n-1; i++)
+            {
+                mean = mean+x[i];
+            }
+            mean = mean/n;
+            
+            //
+            // ADev
+            //
+            for(i=0; i<=n-1; i++)
+            {
+                adev = adev+Math.Abs(x[i]-mean);
+            }
+            adev = adev/n;
+        }
+
+
+        /*************************************************************************
+        Median calculation.
+
+        Input parameters:
+            X   -   sample (array indexes: [0..N-1])
+            N   -   sample size
+
+        Output parameters:
+            Median
+
+          -- ALGLIB --
+             Copyright 06.09.2006 by Bochkanov Sergey
+        *************************************************************************/
+        public static void calculatemedian(double[] x,
+            int n,
+            ref double median)
+        {
+            int i = 0;
+            int ir = 0;
+            int j = 0;
+            int l = 0;
+            int midp = 0;
+            int k = 0;
+            double a = 0;
+            double temp = 0;
+            double tval = 0;
+
+            x = (double[])x.Clone();
+
+            
+            //
+            // Some degenerate cases
+            //
+            median = 0;
+            if( n<=0 )
+            {
+                return;
+            }
+            if( n==1 )
+            {
+                median = x[0];
+                return;
+            }
+            if( n==2 )
+            {
+                median = 0.5*(x[0]+x[1]);
+                return;
+            }
+            
+            //
+            // Common case, N>=3.
+            // Choose X[(N-1)/2]
+            //
+            l = 0;
+            ir = n-1;
+            k = (n-1)/2;
+            while( true )
+            {
+                if( ir<=l+1 )
+                {
+                    
+                    //
+                    // 1 or 2 elements in partition
+                    //
+                    if( ir==l+1 & x[ir]<x[l] )
+                    {
+                        tval = x[l];
+                        x[l] = x[ir];
+                        x[ir] = tval;
+                    }
+                    break;
+                }
+                else
+                {
+                    midp = (l+ir)/2;
+                    tval = x[midp];
+                    x[midp] = x[l+1];
+                    x[l+1] = tval;
+                    if( x[l]>x[ir] )
+                    {
+                        tval = x[l];
+                        x[l] = x[ir];
+                        x[ir] = tval;
+                    }
+                    if( x[l+1]>x[ir] )
+                    {
+                        tval = x[l+1];
+                        x[l+1] = x[ir];
+                        x[ir] = tval;
+                    }
+                    if( x[l]>x[l+1] )
+                    {
+                        tval = x[l];
+                        x[l] = x[l+1];
+                        x[l+1] = tval;
+                    }
+                    i = l+1;
+                    j = ir;
+                    a = x[l+1];
+                    while( true )
+                    {
+                        do
+                        {
+                            i = i+1;
+                        }
+                        while( x[i]<a );
+                        do
+                        {
+                            j = j-1;
+                        }
+                        while( x[j]>a );
+                        if( j<i )
+                        {
+                            break;
+                        }
+                        tval = x[i];
+                        x[i] = x[j];
+                        x[j] = tval;
+                    }
+                    x[l+1] = x[j];
+                    x[j] = a;
+                    if( j>=k )
+                    {
+                        ir = j-1;
+                    }
+                    if( j<=k )
+                    {
+                        l = i;
+                    }
+                }
+            }
+            
+            //
+            // If N is odd, return result
+            //
+            if( n%2==1 )
+            {
+                median = x[k];
+                return;
+            }
+            a = x[n-1];
+            for(i=k+1; i<=n-1; i++)
+            {
+                if( x[i]<a )
+                {
+                    a = x[i];
+                }
+            }
+            median = 0.5*(x[k]+a);
+        }
+
+
+        /*************************************************************************
+        Percentile calculation.
+
+        Input parameters:
+            X   -   sample (array indexes: [0..N-1])
+            N   -   sample size, N>1
+            P   -   percentile (0<=P<=1)
+
+        Output parameters:
+            V   -   percentile
+
+          -- ALGLIB --
+             Copyright 01.03.2008 by Bochkanov Sergey
+        *************************************************************************/
+        public static void calculatepercentile(double[] x,
+            int n,
+            double p,
+            ref double v)
+        {
+            int i1 = 0;
+            double t = 0;
+
+            x = (double[])x.Clone();
+
+            System.Diagnostics.Debug.Assert(n>1, "CalculatePercentile: N<=1!");
+            System.Diagnostics.Debug.Assert(p>=0 & p<=1, "CalculatePercentile: incorrect P!");
+            internalstatheapsort(ref x, n);
+            if( p==0 )
+            {
+                v = x[0];
+                return;
+            }
+            if( p==1 )
+            {
+                v = x[n-1];
+                return;
+            }
+            t = p*(n-1);
+            i1 = (int)Math.Floor(t);
+            t = t-(int)Math.Floor(t);
+            v = x[i1]*(1-t)+x[i1+1]*t;
+        }
+
+
+        private static void internalstatheapsort(ref double[] arr,
+            int n)
+        {
+            int i = 0;
+            int j = 0;
+            int k = 0;
+            int t = 0;
+            double tmp = 0;
+
+            if( n==1 )
+            {
+                return;
+            }
+            i = 2;
+            do
+            {
+                t = i;
+                while( t!=1 )
+                {
+                    k = t/2;
+                    if( arr[k-1]>=arr[t-1] )
+                    {
+                        t = 1;
+                    }
+                    else
+                    {
+                        tmp = arr[k-1];
+                        arr[k-1] = arr[t-1];
+                        arr[t-1] = tmp;
+                        t = k;
+                    }
+                }
+                i = i+1;
+            }
+            while( i<=n );
+            i = n-1;
+            do
+            {
+                tmp = arr[i];
+                arr[i] = arr[0];
+                arr[0] = tmp;
+                t = 1;
+                while( t!=0 )
+                {
+                    k = 2*t;
+                    if( k>i )
+                    {
+                        t = 0;
+                    }
+                    else
+                    {
+                        if( k<i )
+                        {
+                            if( arr[k]>arr[k-1] )
+                            {
+                                k = k+1;
+                            }
+                        }
+                        if( arr[t-1]>=arr[k-1] )
+                        {
+                            t = 0;
+                        }
+                        else
+                        {
+                            tmp = arr[k-1];
+                            arr[k-1] = arr[t-1];
+                            arr[t-1] = tmp;
+                            t = k;
+                        }
+                    }
+                }
+                i = i-1;
+            }
+            while( i>=1 );
+        }
+    }
+}
Index: /trunk/sources/ALGLIB/gammaf.cs
===================================================================
--- /trunk/sources/ALGLIB/gammaf.cs	(revision 2445)
+++ /trunk/sources/ALGLIB/gammaf.cs	(revision 2445)
@@ -0,0 +1,327 @@
+/*************************************************************************
+Cephes Math Library Release 2.8:  June, 2000
+Copyright by Stephen L. Moshier
+
+Contributors:
+    * Sergey Bochkanov (ALGLIB project). Translation from C to
+      pseudocode.
+
+See subroutines comments for additional copyrights.
+
+>>> SOURCE LICENSE >>>
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation (www.fsf.org); either version 2 of the 
+License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+A copy of the GNU General Public License is available at
+http://www.fsf.org/licensing/licenses
+
+>>> END OF LICENSE >>>
+*************************************************************************/
+
+using System;
+
+namespace alglib
+{
+    public class gammaf
+    {
+        /*************************************************************************
+        Gamma function
+
+        Input parameters:
+            X   -   argument
+
+        Domain:
+            0 < X < 171.6
+            -170 < X < 0, X is not an integer.
+
+        Relative error:
+         arithmetic   domain     # trials      peak         rms
+            IEEE    -170,-33      20000       2.3e-15     3.3e-16
+            IEEE     -33,  33     20000       9.4e-16     2.2e-16
+            IEEE      33, 171.6   20000       2.3e-15     3.2e-16
+
+        Cephes Math Library Release 2.8:  June, 2000
+        Original copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
+        Translated to AlgoPascal by Bochkanov Sergey (2005, 2006, 2007).
+        *************************************************************************/
+        public static double gamma(double x)
+        {
+            double result = 0;
+            double p = 0;
+            double pp = 0;
+            double q = 0;
+            double qq = 0;
+            double z = 0;
+            int i = 0;
+            double sgngam = 0;
+
+            sgngam = 1;
+            q = Math.Abs(x);
+            if( q>33.0 )
+            {
+                if( x<0.0 )
+                {
+                    p = (int)Math.Floor(q);
+                    i = (int)Math.Round(p);
+                    if( i%2==0 )
+                    {
+                        sgngam = -1;
+                    }
+                    z = q-p;
+                    if( z>0.5 )
+                    {
+                        p = p+1;
+                        z = q-p;
+                    }
+                    z = q*Math.Sin(Math.PI*z);
+                    z = Math.Abs(z);
+                    z = Math.PI/(z*gammastirf(q));
+                }
+                else
+                {
+                    z = gammastirf(x);
+                }
+                result = sgngam*z;
+                return result;
+            }
+            z = 1;
+            while( x>=3 )
+            {
+                x = x-1;
+                z = z*x;
+            }
+            while( x<0 )
+            {
+                if( x>-0.000000001 )
+                {
+                    result = z/((1+0.5772156649015329*x)*x);
+                    return result;
+                }
+                z = z/x;
+                x = x+1;
+            }
+            while( x<2 )
+            {
+                if( x<0.000000001 )
+                {
+                    result = z/((1+0.5772156649015329*x)*x);
+                    return result;
+                }
+                z = z/x;
+                x = x+1.0;
+            }
+            if( x==2 )
+            {
+                result = z;
+                return result;
+            }
+            x = x-2.0;
+            pp = 1.60119522476751861407E-4;
+            pp = 1.19135147006586384913E-3+x*pp;
+            pp = 1.04213797561761569935E-2+x*pp;
+            pp = 4.76367800457137231464E-2+x*pp;
+            pp = 2.07448227648435975150E-1+x*pp;
+            pp = 4.94214826801497100753E-1+x*pp;
+            pp = 9.99999999999999996796E-1+x*pp;
+            qq = -2.31581873324120129819E-5;
+            qq = 5.39605580493303397842E-4+x*qq;
+            qq = -4.45641913851797240494E-3+x*qq;
+            qq = 1.18139785222060435552E-2+x*qq;
+            qq = 3.58236398605498653373E-2+x*qq;
+            qq = -2.34591795718243348568E-1+x*qq;
+            qq = 7.14304917030273074085E-2+x*qq;
+            qq = 1.00000000000000000320+x*qq;
+            result = z*pp/qq;
+            return result;
+            return result;
+        }
+
+
+        /*************************************************************************
+        Natural logarithm of gamma function
+
+        Input parameters:
+            X       -   argument
+
+        Result:
+            logarithm of the absolute value of the Gamma(X).
+
+        Output parameters:
+            SgnGam  -   sign(Gamma(X))
+
+        Domain:
+            0 < X < 2.55e305
+            -2.55e305 < X < 0, X is not an integer.
+
+        ACCURACY:
+        arithmetic      domain        # trials     peak         rms
+           IEEE    0, 3                 28000     5.4e-16     1.1e-16
+           IEEE    2.718, 2.556e305     40000     3.5e-16     8.3e-17
+        The error criterion was relative when the function magnitude
+        was greater than one but absolute when it was less than one.
+
+        The following test used the relative error criterion, though
+        at certain points the relative error could be much higher than
+        indicated.
+           IEEE    -200, -4             10000     4.8e-16     1.3e-16
+
+        Cephes Math Library Release 2.8:  June, 2000
+        Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
+        Translated to AlgoPascal by Bochkanov Sergey (2005, 2006, 2007).
+        *************************************************************************/
+        public static double lngamma(double x,
+            ref double sgngam)
+        {
+            double result = 0;
+            double a = 0;
+            double b = 0;
+            double c = 0;
+            double p = 0;
+            double q = 0;
+            double u = 0;
+            double w = 0;
+            double z = 0;
+            int i = 0;
+            double logpi = 0;
+            double ls2pi = 0;
+            double tmp = 0;
+
+            sgngam = 1;
+            logpi = 1.14472988584940017414;
+            ls2pi = 0.91893853320467274178;
+            if( x<-34.0 )
+            {
+                q = -x;
+                w = lngamma(q, ref tmp);
+                p = (int)Math.Floor(q);
+                i = (int)Math.Round(p);
+                if( i%2==0 )
+                {
+                    sgngam = -1;
+                }
+                else
+                {
+                    sgngam = 1;
+                }
+                z = q-p;
+                if( z>0.5 )
+                {
+                    p = p+1;
+                    z = p-q;
+                }
+                z = q*Math.Sin(Math.PI*z);
+                result = logpi-Math.Log(z)-w;
+                return result;
+            }
+            if( x<13 )
+            {
+                z = 1;
+                p = 0;
+                u = x;
+                while( u>=3 )
+                {
+                    p = p-1;
+                    u = x+p;
+                    z = z*u;
+                }
+                while( u<2 )
+                {
+                    z = z/u;
+                    p = p+1;
+                    u = x+p;
+                }
+                if( z<0 )
+                {
+                    sgngam = -1;
+                    z = -z;
+                }
+                else
+                {
+                    sgngam = 1;
+                }
+                if( u==2 )
+                {
+                    result = Math.Log(z);
+                    return result;
+                }
+                p = p-2;
+                x = x+p;
+                b = -1378.25152569120859100;
+                b = -38801.6315134637840924+x*b;
+                b = -331612.992738871184744+x*b;
+                b = -1162370.97492762307383+x*b;
+                b = -1721737.00820839662146+x*b;
+                b = -853555.664245765465627+x*b;
+                c = 1;
+                c = -351.815701436523470549+x*c;
+                c = -17064.2106651881159223+x*c;
+                c = -220528.590553854454839+x*c;
+                c = -1139334.44367982507207+x*c;
+                c = -2532523.07177582951285+x*c;
+                c = -2018891.41433532773231+x*c;
+                p = x*b/c;
+                result = Math.Log(z)+p;
+                return result;
+            }
+            q = (x-0.5)*Math.Log(x)-x+ls2pi;
+            if( x>100000000 )
+            {
+                result = q;
+                return result;
+            }
+            p = 1/(x*x);
+            if( x>=1000.0 )
+            {
+                q = q+((7.9365079365079365079365*0.0001*p-2.7777777777777777777778*0.001)*p+0.0833333333333333333333)/x;
+            }
+            else
+            {
+                a = 8.11614167470508450300*0.0001;
+                a = -(5.95061904284301438324*0.0001)+p*a;
+                a = 7.93650340457716943945*0.0001+p*a;
+                a = -(2.77777777730099687205*0.001)+p*a;
+                a = 8.33333333333331927722*0.01+p*a;
+                q = q+a/x;
+            }
+            result = q;
+            return result;
+        }
+
+
+        private static double gammastirf(double x)
+        {
+            double result = 0;
+            double y = 0;
+            double w = 0;
+            double v = 0;
+            double stir = 0;
+
+            w = 1/x;
+            stir = 7.87311395793093628397E-4;
+            stir = -2.29549961613378126380E-4+w*stir;
+            stir = -2.68132617805781232825E-3+w*stir;
+            stir = 3.47222221605458667310E-3+w*stir;
+            stir = 8.33333333333482257126E-2+w*stir;
+            w = 1+w*stir;
+            y = Math.Exp(x);
+            if( x>143.01608 )
+            {
+                v = Math.Pow(x, 0.5*x-0.25);
+                y = v*(v/y);
+            }
+            else
+            {
+                y = Math.Pow(x, x-0.5)/y;
+            }
+            result = 2.50662827463100050242*y*w;
+            return result;
+        }
+    }
+}
Index: /trunk/sources/ALGLIB/igammaf.cs
===================================================================
--- /trunk/sources/ALGLIB/igammaf.cs	(revision 2445)
+++ /trunk/sources/ALGLIB/igammaf.cs	(revision 2445)
@@ -0,0 +1,423 @@
+/*************************************************************************
+Cephes Math Library Release 2.8:  June, 2000
+Copyright by Stephen L. Moshier
+
+Contributors:
+    * Sergey Bochkanov (ALGLIB project). Translation from C to
+      pseudocode.
+
+See subroutines comments for additional copyrights.
+
+>>> SOURCE LICENSE >>>
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation (www.fsf.org); either version 2 of the 
+License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+A copy of the GNU General Public License is available at
+http://www.fsf.org/licensing/licenses
+
+>>> END OF LICENSE >>>
+*************************************************************************/
+
+using System;
+
+namespace alglib
+{
+    public class igammaf
+    {
+        /*************************************************************************
+        Incomplete gamma integral
+
+        The function is defined by
+
+                                  x
+                                   -
+                          1       | |  -t  a-1
+         igam(a,x)  =   -----     |   e   t   dt.
+                         -      | |
+                        | (a)    -
+                                  0
+
+
+        In this implementation both arguments must be positive.
+        The integral is evaluated by either a power series or
+        continued fraction expansion, depending on the relative
+        values of a and x.
+
+        ACCURACY:
+
+                             Relative error:
+        arithmetic   domain     # trials      peak         rms
+           IEEE      0,30       200000       3.6e-14     2.9e-15
+           IEEE      0,100      300000       9.9e-14     1.5e-14
+
+        Cephes Math Library Release 2.8:  June, 2000
+        Copyright 1985, 1987, 2000 by Stephen L. Moshier
+        *************************************************************************/
+        public static double incompletegamma(double a,
+            double x)
+        {
+            double result = 0;
+            double igammaepsilon = 0;
+            double ans = 0;
+            double ax = 0;
+            double c = 0;
+            double r = 0;
+            double tmp = 0;
+
+            igammaepsilon = 0.000000000000001;
+            if( x<=0 | a<=0 )
+            {
+                result = 0;
+                return result;
+            }
+            if( x>1 & x>a )
+            {
+                result = 1-incompletegammac(a, x);
+                return result;
+            }
+            ax = a*Math.Log(x)-x-gammaf.lngamma(a, ref tmp);
+            if( ax<-709.78271289338399 )
+            {
+                result = 0;
+                return result;
+            }
+            ax = Math.Exp(ax);
+            r = a;
+            c = 1;
+            ans = 1;
+            do
+            {
+                r = r+1;
+                c = c*x/r;
+                ans = ans+c;
+            }
+            while( c/ans>igammaepsilon );
+            result = ans*ax/a;
+            return result;
+        }
+
+
+        /*************************************************************************
+        Complemented incomplete gamma integral
+
+        The function is defined by
+
+
+         igamc(a,x)   =   1 - igam(a,x)
+
+                                   inf.
+                                     -
+                            1       | |  -t  a-1
+                      =   -----     |   e   t   dt.
+                           -      | |
+                          | (a)    -
+                                    x
+
+
+        In this implementation both arguments must be positive.
+        The integral is evaluated by either a power series or
+        continued fraction expansion, depending on the relative
+        values of a and x.
+
+        ACCURACY:
+
+        Tested at random a, x.
+                       a         x                      Relative error:
+        arithmetic   domain   domain     # trials      peak         rms
+           IEEE     0.5,100   0,100      200000       1.9e-14     1.7e-15
+           IEEE     0.01,0.5  0,100      200000       1.4e-13     1.6e-15
+
+        Cephes Math Library Release 2.8:  June, 2000
+        Copyright 1985, 1987, 2000 by Stephen L. Moshier
+        *************************************************************************/
+        public static double incompletegammac(double a,
+            double x)
+        {
+            double result = 0;
+            double igammaepsilon = 0;
+            double igammabignumber = 0;
+            double igammabignumberinv = 0;
+            double ans = 0;
+            double ax = 0;
+            double c = 0;
+            double yc = 0;
+            double r = 0;
+            double t = 0;
+            double y = 0;
+            double z = 0;
+            double pk = 0;
+            double pkm1 = 0;
+            double pkm2 = 0;
+            double qk = 0;
+            double qkm1 = 0;
+            double qkm2 = 0;
+            double tmp = 0;
+
+            igammaepsilon = 0.000000000000001;
+            igammabignumber = 4503599627370496.0;
+            igammabignumberinv = 2.22044604925031308085*0.0000000000000001;
+            if( x<=0 | a<=0 )
+            {
+                result = 1;
+                return result;
+            }
+            if( x<1 | x<a )
+            {
+                result = 1-incompletegamma(a, x);
+                return result;
+            }
+            ax = a*Math.Log(x)-x-gammaf.lngamma(a, ref tmp);
+            if( ax<-709.78271289338399 )
+            {
+                result = 0;
+                return result;
+            }
+            ax = Math.Exp(ax);
+            y = 1-a;
+            z = x+y+1;
+            c = 0;
+            pkm2 = 1;
+            qkm2 = x;
+            pkm1 = x+1;
+            qkm1 = z*x;
+            ans = pkm1/qkm1;
+            do
+            {
+                c = c+1;
+                y = y+1;
+                z = z+2;
+                yc = y*c;
+                pk = pkm1*z-pkm2*yc;
+                qk = qkm1*z-qkm2*yc;
+                if( qk!=0 )
+                {
+                    r = pk/qk;
+                    t = Math.Abs((ans-r)/r);
+                    ans = r;
+                }
+                else
+                {
+                    t = 1;
+                }
+                pkm2 = pkm1;
+                pkm1 = pk;
+                qkm2 = qkm1;
+                qkm1 = qk;
+                if( Math.Abs(pk)>igammabignumber )
+                {
+                    pkm2 = pkm2*igammabignumberinv;
+                    pkm1 = pkm1*igammabignumberinv;
+                    qkm2 = qkm2*igammabignumberinv;
+                    qkm1 = qkm1*igammabignumberinv;
+                }
+            }
+            while( t>igammaepsilon );
+            result = ans*ax;
+            return result;
+        }
+
+
+        /*************************************************************************
+        Inverse of complemented imcomplete gamma integral
+
+        Given p, the function finds x such that
+
+         igamc( a, x ) = p.
+
+        Starting with the approximate value
+
+                3
+         x = a t
+
+         where
+
+         t = 1 - d - ndtri(p) sqrt(d)
+
+        and
+
+         d = 1/9a,
+
+        the routine performs up to 10 Newton iterations to find the
+        root of igamc(a,x) - p = 0.
+
+        ACCURACY:
+
+        Tested at random a, p in the intervals indicated.
+
+                       a        p                      Relative error:
+        arithmetic   domain   domain     # trials      peak         rms
+           IEEE     0.5,100   0,0.5       100000       1.0e-14     1.7e-15
+           IEEE     0.01,0.5  0,0.5       100000       9.0e-14     3.4e-15
+           IEEE    0.5,10000  0,0.5        20000       2.3e-13     3.8e-14
+
+        Cephes Math Library Release 2.8:  June, 2000
+        Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
+        *************************************************************************/
+        public static double invincompletegammac(double a,
+            double y0)
+        {
+            double result = 0;
+            double igammaepsilon = 0;
+            double iinvgammabignumber = 0;
+            double x0 = 0;
+            double x1 = 0;
+            double x = 0;
+            double yl = 0;
+            double yh = 0;
+            double y = 0;
+            double d = 0;
+            double lgm = 0;
+            double dithresh = 0;
+            int i = 0;
+            int dir = 0;
+            double tmp = 0;
+
+            igammaepsilon = 0.000000000000001;
+            iinvgammabignumber = 4503599627370496.0;
+            x0 = iinvgammabignumber;
+            yl = 0;
+            x1 = 0;
+            yh = 1;
+            dithresh = 5*igammaepsilon;
+            d = 1/(9*a);
+            y = 1-d-normaldistr.invnormaldistribution(y0)*Math.Sqrt(d);
+            x = a*y*y*y;
+            lgm = gammaf.lngamma(a, ref tmp);
+            i = 0;
+            while( i<10 )
+            {
+                if( x>x0 | x<x1 )
+                {
+                    d = 0.0625;
+                    break;
+                }
+                y = incompletegammac(a, x);
+                if( y<yl | y>yh )
+                {
+                    d = 0.0625;
+                    break;
+                }
+                if( y<y0 )
+                {
+                    x0 = x;
+                    yl = y;
+                }
+                else
+                {
+                    x1 = x;
+                    yh = y;
+                }
+                d = (a-1)*Math.Log(x)-x-lgm;
+                if( d<-709.78271289338399 )
+                {
+                    d = 0.0625;
+                    break;
+                }
+                d = -Math.Exp(d);
+                d = (y-y0)/d;
+                if( Math.Abs(d/x)<igammaepsilon )
+                {
+                    result = x;
+                    return result;
+                }
+                x = x-d;
+                i = i+1;
+            }
+            if( x0==iinvgammabignumber )
+            {
+                if( x<=0 )
+                {
+                    x = 1;
+                }
+                while( x0==iinvgammabignumber )
+                {
+                    x = (1+d)*x;
+                    y = incompletegammac(a, x);
+                    if( y<y0 )
+                    {
+                        x0 = x;
+                        yl = y;
+                        break;
+                    }
+                    d = d+d;
+                }
+            }
+            d = 0.5;
+            dir = 0;
+            i = 0;
+            while( i<400 )
+            {
+                x = x1+d*(x0-x1);
+                y = incompletegammac(a, x);
+                lgm = (x0-x1)/(x1+x0);
+                if( Math.Abs(lgm)<dithresh )
+                {
+                    break;
+                }
+                lgm = (y-y0)/y0;
+                if( Math.Abs(lgm)<dithresh )
+                {
+                    break;
+                }
+                if( x<=0.0 )
+                {
+                    break;
+                }
+                if( y>=y0 )
+                {
+                    x1 = x;
+                    yh = y;
+                    if( dir<0 )
+                    {
+                        dir = 0;
+                        d = 0.5;
+                    }
+                    else
+                    {
+                        if( dir>1 )
+                        {
+                            d = 0.5*d+0.5;
+                        }
+                        else
+                        {
+                            d = (y0-yl)/(yh-yl);
+                        }
+                    }
+                    dir = dir+1;
+                }
+                else
+                {
+                    x0 = x;
+                    yl = y;
+                    if( dir>0 )
+                    {
+                        dir = 0;
+                        d = 0.5;
+                    }
+                    else
+                    {
+                        if( dir<-1 )
+                        {
+                            d = 0.5*d;
+                        }
+                        else
+                        {
+                            d = (y0-yl)/(yh-yl);
+                        }
+                    }
+                    dir = dir-1;
+                }
+                i = i+1;
+            }
+            result = x;
+            return result;
+        }
+    }
+}
Index: /trunk/sources/ALGLIB/linreg.cs
===================================================================
--- /trunk/sources/ALGLIB/linreg.cs	(revision 2445)
+++ /trunk/sources/ALGLIB/linreg.cs	(revision 2445)
@@ -0,0 +1,1348 @@
+/*************************************************************************
+Copyright (c) 2007-2008, Sergey Bochkanov (ALGLIB project).
+
+>>> SOURCE LICENSE >>>
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation (www.fsf.org); either version 2 of the 
+License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+A copy of the GNU General Public License is available at
+http://www.fsf.org/licensing/licenses
+
+>>> END OF LICENSE >>>
+*************************************************************************/
+
+using System;
+
+namespace alglib
+{
+    public class linreg
+    {
+        public struct linearmodel
+        {
+            public double[] w;
+        };
+
+
+        /*************************************************************************
+        LRReport structure contains additional information about linear model:
+        * C             -   covariation matrix,  array[0..NVars,0..NVars].
+                            C[i,j] = Cov(A[i],A[j])
+        * RMSError      -   root mean square error on a training set
+        * AvgError      -   average error on a training set
+        * AvgRelError   -   average relative error on a training set (excluding
+                            observations with zero function value).
+        * CVRMSError    -   leave-one-out cross-validation estimate of
+                            generalization error. Calculated using fast algorithm
+                            with O(NVars*NPoints) complexity.
+        * CVAvgError    -   cross-validation estimate of average error
+        * CVAvgRelError -   cross-validation estimate of average relative error
+
+        All other fields of the structure are intended for internal use and should
+        not be used outside ALGLIB.
+        *************************************************************************/
+        public struct lrreport
+        {
+            public double[,] c;
+            public double rmserror;
+            public double avgerror;
+            public double avgrelerror;
+            public double cvrmserror;
+            public double cvavgerror;
+            public double cvavgrelerror;
+            public int ncvdefects;
+            public int[] cvdefects;
+        };
+
+
+
+
+        public const int lrvnum = 5;
+
+
+        /*************************************************************************
+        Linear regression
+
+        Subroutine builds model:
+
+            Y = A(0)*X[0] + ... + A(N-1)*X[N-1] + A(N)
+
+        and model found in ALGLIB format, covariation matrix, training set  errors
+        (rms,  average,  average  relative)   and  leave-one-out  cross-validation
+        estimate of the generalization error. CV  estimate calculated  using  fast
+        algorithm with O(NPoints*NVars) complexity.
+
+        When  covariation  matrix  is  calculated  standard deviations of function
+        values are assumed to be equal to RMS error on the training set.
+
+        INPUT PARAMETERS:
+            XY          -   training set, array [0..NPoints-1,0..NVars]:
+                            * NVars columns - independent variables
+                            * last column - dependent variable
+            NPoints     -   training set size, NPoints>NVars+1
+            NVars       -   number of independent variables
+
+        OUTPUT PARAMETERS:
+            Info        -   return code:
+                            * -255, in case of unknown internal error
+                            * -4, if internal SVD subroutine haven't converged
+                            * -1, if incorrect parameters was passed (NPoints<NVars+2, NVars<1).
+                            *  1, if subroutine successfully finished
+            LM          -   linear model in the ALGLIB format. Use subroutines of
+                            this unit to work with the model.
+            AR          -   additional results
+
+
+          -- ALGLIB --
+             Copyright 02.08.2008 by Bochkanov Sergey
+        *************************************************************************/
+        public static void lrbuild(ref double[,] xy,
+            int npoints,
+            int nvars,
+            ref int info,
+            ref linearmodel lm,
+            ref lrreport ar)
+        {
+            double[] s = new double[0];
+            int i = 0;
+            double sigma2 = 0;
+            int i_ = 0;
+
+            if( npoints<=nvars+1 | nvars<1 )
+            {
+                info = -1;
+                return;
+            }
+            s = new double[npoints-1+1];
+            for(i=0; i<=npoints-1; i++)
+            {
+                s[i] = 1;
+            }
+            lrbuilds(ref xy, ref s, npoints, nvars, ref info, ref lm, ref ar);
+            if( info<0 )
+            {
+                return;
+            }
+            sigma2 = AP.Math.Sqr(ar.rmserror)*npoints/(npoints-nvars-1);
+            for(i=0; i<=nvars; i++)
+            {
+                for(i_=0; i_<=nvars;i_++)
+                {
+                    ar.c[i,i_] = sigma2*ar.c[i,i_];
+                }
+            }
+        }
+
+
+        /*************************************************************************
+        Linear regression
+
+        Variant of LRBuild which uses vector of standatd deviations (errors in
+        function values).
+
+        INPUT PARAMETERS:
+            XY          -   training set, array [0..NPoints-1,0..NVars]:
+                            * NVars columns - independent variables
+                            * last column - dependent variable
+            S           -   standard deviations (errors in function values)
+                            array[0..NPoints-1], S[i]>0.
+            NPoints     -   training set size, NPoints>NVars+1
+            NVars       -   number of independent variables
+
+        OUTPUT PARAMETERS:
+            Info        -   return code:
+                            * -255, in case of unknown internal error
+                            * -4, if internal SVD subroutine haven't converged
+                            * -1, if incorrect parameters was passed (NPoints<NVars+2, NVars<1).
+                            * -2, if S[I]<=0
+                            *  1, if subroutine successfully finished
+            LM          -   linear model in the ALGLIB format. Use subroutines of
+                            this unit to work with the model.
+            AR          -   additional results
+
+
+          -- ALGLIB --
+             Copyright 02.08.2008 by Bochkanov Sergey
+        *************************************************************************/
+        public static void lrbuilds(ref double[,] xy,
+            ref double[] s,
+            int npoints,
+            int nvars,
+            ref int info,
+            ref linearmodel lm,
+            ref lrreport ar)
+        {
+            double[,] xyi = new double[0,0];
+            double[] x = new double[0];
+            double[] means = new double[0];
+            double[] sigmas = new double[0];
+            int i = 0;
+            int j = 0;
+            double v = 0;
+            int offs = 0;
+            double mean = 0;
+            double variance = 0;
+            double skewness = 0;
+            double kurtosis = 0;
+            int i_ = 0;
+
+            
+            //
+            // Test parameters
+            //
+            if( npoints<=nvars+1 | nvars<1 )
+            {
+                info = -1;
+                return;
+            }
+            
+            //
+            // Copy data, add one more column (constant term)
+            //
+            xyi = new double[npoints-1+1, nvars+1+1];
+            for(i=0; i<=npoints-1; i++)
+            {
+                for(i_=0; i_<=nvars-1;i_++)
+                {
+                    xyi[i,i_] = xy[i,i_];
+                }
+                xyi[i,nvars] = 1;
+                xyi[i,nvars+1] = xy[i,nvars];
+            }
+            
+            //
+            // Standartization
+            //
+            x = new double[npoints-1+1];
+            means = new double[nvars-1+1];
+            sigmas = new double[nvars-1+1];
+            for(j=0; j<=nvars-1; j++)
+            {
+                for(i_=0; i_<=npoints-1;i_++)
+                {
+                    x[i_] = xy[i_,j];
+                }
+                descriptivestatistics.calculatemoments(ref x, npoints, ref mean, ref variance, ref skewness, ref kurtosis);
+                means[j] = mean;
+                sigmas[j] = Math.Sqrt(variance);
+                if( sigmas[j]==0 )
+                {
+                    sigmas[j] = 1;
+                }
+                for(i=0; i<=npoints-1; i++)
+                {
+                    xyi[i,j] = (xyi[i,j]-means[j])/sigmas[j];
+                }
+            }
+            
+            //
+            // Internal processing
+            //
+            lrinternal(ref xyi, ref s, npoints, nvars+1, ref info, ref lm, ref ar);
+            if( info<0 )
+            {
+                return;
+            }
+            
+            //
+            // Un-standartization
+            //
+            offs = (int)Math.Round(lm.w[3]);
+            for(j=0; j<=nvars-1; j++)
+            {
+                
+                //
+                // Constant term is updated (and its covariance too,
+                // since it gets some variance from J-th component)
+                //
+                lm.w[offs+nvars] = lm.w[offs+nvars]-lm.w[offs+j]*means[j]/sigmas[j];
+                v = means[j]/sigmas[j];
+                for(i_=0; i_<=nvars;i_++)
+                {
+                    ar.c[nvars,i_] = ar.c[nvars,i_] - v*ar.c[j,i_];
+                }
+                for(i_=0; i_<=nvars;i_++)
+                {
+                    ar.c[i_,nvars] = ar.c[i_,nvars] - v*ar.c[i_,j];
+                }
+                
+                //
+                // J-th term is updated
+                //
+                lm.w[offs+j] = lm.w[offs+j]/sigmas[j];
+                v = 1/sigmas[j];
+                for(i_=0; i_<=nvars;i_++)
+                {
+                    ar.c[j,i_] = v*ar.c[j,i_];
+                }
+                for(i_=0; i_<=nvars;i_++)
+                {
+                    ar.c[i_,j] = v*ar.c[i_,j];
+                }
+            }
+        }
+
+
+        /*************************************************************************
+        Like LRBuildS, but builds model
+
+            Y = A(0)*X[0] + ... + A(N-1)*X[N-1]
+
+        i.e. with zero constant term.
+
+          -- ALGLIB --
+             Copyright 30.10.2008 by Bochkanov Sergey
+        *************************************************************************/
+        public static void lrbuildzs(ref double[,] xy,
+            ref double[] s,
+            int npoints,
+            int nvars,
+            ref int info,
+            ref linearmodel lm,
+            ref lrreport ar)
+        {
+            double[,] xyi = new double[0,0];
+            double[] x = new double[0];
+            double[] c = new double[0];
+            int i = 0;
+            int j = 0;
+            double v = 0;
+            int offs = 0;
+            double mean = 0;
+            double variance = 0;
+            double skewness = 0;
+            double kurtosis = 0;
+            int i_ = 0;
+
+            
+            //
+            // Test parameters
+            //
+            if( npoints<=nvars+1 | nvars<1 )
+            {
+                info = -1;
+                return;
+            }
+            
+            //
+            // Copy data, add one more column (constant term)
+            //
+            xyi = new double[npoints-1+1, nvars+1+1];
+            for(i=0; i<=npoints-1; i++)
+            {
+                for(i_=0; i_<=nvars-1;i_++)
+                {
+                    xyi[i,i_] = xy[i,i_];
+                }
+                xyi[i,nvars] = 0;
+                xyi[i,nvars+1] = xy[i,nvars];
+            }
+            
+            //
+            // Standartization: unusual scaling
+            //
+            x = new double[npoints-1+1];
+            c = new double[nvars-1+1];
+            for(j=0; j<=nvars-1; j++)
+            {
+                for(i_=0; i_<=npoints-1;i_++)
+                {
+                    x[i_] = xy[i_,j];
+                }
+                descriptivestatistics.calculatemoments(ref x, npoints, ref mean, ref variance, ref skewness, ref kurtosis);
+                if( Math.Abs(mean)>Math.Sqrt(variance) )
+                {
+                    
+                    //
+                    // variation is relatively small, it is better to
+                    // bring mean value to 1
+                    //
+                    c[j] = mean;
+                }
+                else
+                {
+                    
+                    //
+                    // variation is large, it is better to bring variance to 1
+                    //
+                    if( variance==0 )
+                    {
+                        variance = 1;
+                    }
+                    c[j] = Math.Sqrt(variance);
+                }
+                for(i=0; i<=npoints-1; i++)
+                {
+                    xyi[i,j] = xyi[i,j]/c[j];
+                }
+            }
+            
+            //
+            // Internal processing
+            //
+            lrinternal(ref xyi, ref s, npoints, nvars+1, ref info, ref lm, ref ar);
+            if( info<0 )
+            {
+                return;
+            }
+            
+            //
+            // Un-standartization
+            //
+            offs = (int)Math.Round(lm.w[3]);
+            for(j=0; j<=nvars-1; j++)
+            {
+                
+                //
+                // J-th term is updated
+                //
+                lm.w[offs+j] = lm.w[offs+j]/c[j];
+                v = 1/c[j];
+                for(i_=0; i_<=nvars;i_++)
+                {
+                    ar.c[j,i_] = v*ar.c[j,i_];
+                }
+                for(i_=0; i_<=nvars;i_++)
+                {
+                    ar.c[i_,j] = v*ar.c[i_,j];
+                }
+            }
+        }
+
+
+        /*************************************************************************
+        Like LRBuild but builds model
+
+            Y = A(0)*X[0] + ... + A(N-1)*X[N-1]
+
+        i.e. with zero constant term.
+
+          -- ALGLIB --
+             Copyright 30.10.2008 by Bochkanov Sergey
+        *************************************************************************/
+        public static void lrbuildz(ref double[,] xy,
+            int npoints,
+            int nvars,
+            ref int info,
+            ref linearmodel lm,
+            ref lrreport ar)
+        {
+            double[] s = new double[0];
+            int i = 0;
+            double sigma2 = 0;
+            int i_ = 0;
+
+            if( npoints<=nvars+1 | nvars<1 )
+            {
+                info = -1;
+                return;
+            }
+            s = new double[npoints-1+1];
+            for(i=0; i<=npoints-1; i++)
+            {
+                s[i] = 1;
+            }
+            lrbuildzs(ref xy, ref s, npoints, nvars, ref info, ref lm, ref ar);
+            if( info<0 )
+            {
+                return;
+            }
+            sigma2 = AP.Math.Sqr(ar.rmserror)*npoints/(npoints-nvars-1);
+            for(i=0; i<=nvars; i++)
+            {
+                for(i_=0; i_<=nvars;i_++)
+                {
+                    ar.c[i,i_] = sigma2*ar.c[i,i_];
+                }
+            }
+        }
+
+
+        /*************************************************************************
+        Unpacks coefficients of linear model.
+
+        INPUT PARAMETERS:
+            LM          -   linear model in ALGLIB format
+
+        OUTPUT PARAMETERS:
+            V           -   coefficients, array[0..NVars]
+            NVars       -   number of independent variables (one less than number
+                            of coefficients)
+
+          -- ALGLIB --
+             Copyright 30.08.2008 by Bochkanov Sergey
+        *************************************************************************/
+        public static void lrunpack(ref linearmodel lm,
+            ref double[] v,
+            ref int nvars)
+        {
+            int offs = 0;
+            int i_ = 0;
+            int i1_ = 0;
+
+            System.Diagnostics.Debug.Assert((int)Math.Round(lm.w[1])==lrvnum, "LINREG: Incorrect LINREG version!");
+            nvars = (int)Math.Round(lm.w[2]);
+            offs = (int)Math.Round(lm.w[3]);
+            v = new double[nvars+1];
+            i1_ = (offs) - (0);
+            for(i_=0; i_<=nvars;i_++)
+            {
+                v[i_] = lm.w[i_+i1_];
+            }
+        }
+
+
+        /*************************************************************************
+        "Packs" coefficients and creates linear model in ALGLIB format (LRUnpack
+        reversed).
+
+        INPUT PARAMETERS:
+            V           -   coefficients, array[0..NVars]
+            NVars       -   number of independent variables
+
+        OUTPUT PAREMETERS:
+            LM          -   linear model.
+
+          -- ALGLIB --
+             Copyright 30.08.2008 by Bochkanov Sergey
+        *************************************************************************/
+        public static void lrpack(ref double[] v,
+            int nvars,
+            ref linearmodel lm)
+        {
+            int offs = 0;
+            int i_ = 0;
+            int i1_ = 0;
+
+            lm.w = new double[4+nvars+1];
+            offs = 4;
+            lm.w[0] = 4+nvars+1;
+            lm.w[1] = lrvnum;
+            lm.w[2] = nvars;
+            lm.w[3] = offs;
+            i1_ = (0) - (offs);
+            for(i_=offs; i_<=offs+nvars;i_++)
+            {
+                lm.w[i_] = v[i_+i1_];
+            }
+        }
+
+
+        /*************************************************************************
+        Procesing
+
+        INPUT PARAMETERS:
+            LM      -   linear model
+            X       -   input vector,  array[0..NVars-1].
+
+        Result:
+            value of linear model regression estimate
+
+          -- ALGLIB --
+             Copyright 03.09.2008 by Bochkanov Sergey
+        *************************************************************************/
+        public static double lrprocess(ref linearmodel lm,
+            ref double[] x)
+        {
+            double result = 0;
+            double v = 0;
+            int offs = 0;
+            int nvars = 0;
+            int i_ = 0;
+            int i1_ = 0;
+
+            System.Diagnostics.Debug.Assert((int)Math.Round(lm.w[1])==lrvnum, "LINREG: Incorrect LINREG version!");
+            nvars = (int)Math.Round(lm.w[2]);
+            offs = (int)Math.Round(lm.w[3]);
+            i1_ = (offs)-(0);
+            v = 0.0;
+            for(i_=0; i_<=nvars-1;i_++)
+            {
+                v += x[i_]*lm.w[i_+i1_];
+            }
+            result = v+lm.w[offs+nvars];
+            return result;
+        }
+
+
+        /*************************************************************************
+        RMS error on the test set
+
+        INPUT PARAMETERS:
+            LM      -   linear model
+            XY      -   test set
+            NPoints -   test set size
+
+        RESULT:
+            root mean square error.
+
+          -- ALGLIB --
+             Copyright 30.08.2008 by Bochkanov Sergey
+        *************************************************************************/
+        public static double lrrmserror(ref linearmodel lm,
+            ref double[,] xy,
+            int npoints)
+        {
+            double result = 0;
+            int i = 0;
+            double v = 0;
+            int offs = 0;
+            int nvars = 0;
+            int i_ = 0;
+            int i1_ = 0;
+
+            System.Diagnostics.Debug.Assert((int)Math.Round(lm.w[1])==lrvnum, "LINREG: Incorrect LINREG version!");
+            nvars = (int)Math.Round(lm.w[2]);
+            offs = (int)Math.Round(lm.w[3]);
+            result = 0;
+            for(i=0; i<=npoints-1; i++)
+            {
+                i1_ = (offs)-(0);
+                v = 0.0;
+                for(i_=0; i_<=nvars-1;i_++)
+                {
+                    v += xy[i,i_]*lm.w[i_+i1_];
+                }
+                v = v+lm.w[offs+nvars];
+                result = result+AP.Math.Sqr(v-xy[i,nvars]);
+            }
+            result = Math.Sqrt(result/npoints);
+            return result;
+        }
+
+
+        /*************************************************************************
+        Average error on the test set
+
+        INPUT PARAMETERS:
+            LM      -   linear model
+            XY      -   test set
+            NPoints -   test set size
+
+        RESULT:
+            average error.
+
+          -- ALGLIB --
+             Copyright 30.08.2008 by Bochkanov Sergey
+        *************************************************************************/
+        public static double lravgerror(ref linearmodel lm,
+            ref double[,] xy,
+            int npoints)
+        {
+            double result = 0;
+            int i = 0;
+            double v = 0;
+            int offs = 0;
+            int nvars = 0;
+            int i_ = 0;
+            int i1_ = 0;
+
+            System.Diagnostics.Debug.Assert((int)Math.Round(lm.w[1])==lrvnum, "LINREG: Incorrect LINREG version!");
+            nvars = (int)Math.Round(lm.w[2]);
+            offs = (int)Math.Round(lm.w[3]);
+            result = 0;
+            for(i=0; i<=npoints-1; i++)
+            {
+                i1_ = (offs)-(0);
+                v = 0.0;
+                for(i_=0; i_<=nvars-1;i_++)
+                {
+                    v += xy[i,i_]*lm.w[i_+i1_];
+                }
+                v = v+lm.w[offs+nvars];
+                result = result+Math.Abs(v-xy[i,nvars]);
+            }
+            result = result/npoints;
+            return result;
+        }
+
+
+        /*************************************************************************
+        RMS error on the test set
+
+        INPUT PARAMETERS:
+            LM      -   linear model
+            XY      -   test set
+            NPoints -   test set size
+
+        RESULT:
+            average relative error.
+
+          -- ALGLIB --
+             Copyright 30.08.2008 by Bochkanov Sergey
+        *************************************************************************/
+        public static double lravgrelerror(ref linearmodel lm,
+            ref double[,] xy,
+            int npoints)
+        {
+            double result = 0;
+            int i = 0;
+            int k = 0;
+            double v = 0;
+            int offs = 0;
+            int nvars = 0;
+            int i_ = 0;
+            int i1_ = 0;
+
+            System.Diagnostics.Debug.Assert((int)Math.Round(lm.w[1])==lrvnum, "LINREG: Incorrect LINREG version!");
+            nvars = (int)Math.Round(lm.w[2]);
+            offs = (int)Math.Round(lm.w[3]);
+            result = 0;
+            k = 0;
+            for(i=0; i<=npoints-1; i++)
+            {
+                if( xy[i,nvars]!=0 )
+                {
+                    i1_ = (offs)-(0);
+                    v = 0.0;
+                    for(i_=0; i_<=nvars-1;i_++)
+                    {
+                        v += xy[i,i_]*lm.w[i_+i1_];
+                    }
+                    v = v+lm.w[offs+nvars];
+                    result = result+Math.Abs((v-xy[i,nvars])/xy[i,nvars]);
+                    k = k+1;
+                }
+            }
+            if( k!=0 )
+            {
+                result = result/k;
+            }
+            return result;
+        }
+
+
+        /*************************************************************************
+        Copying of LinearModel strucure
+
+        INPUT PARAMETERS:
+            LM1 -   original
+
+        OUTPUT PARAMETERS:
+            LM2 -   copy
+
+          -- ALGLIB --
+             Copyright 15.03.2009 by Bochkanov Sergey
+        *************************************************************************/
+        public static void lrcopy(ref linearmodel lm1,
+            ref linearmodel lm2)
+        {
+            int k = 0;
+            int i_ = 0;
+
+            k = (int)Math.Round(lm1.w[0]);
+            lm2.w = new double[k-1+1];
+            for(i_=0; i_<=k-1;i_++)
+            {
+                lm2.w[i_] = lm1.w[i_];
+            }
+        }
+
+
+        /*************************************************************************
+        Serialization of LinearModel strucure
+
+        INPUT PARAMETERS:
+            LM      -   original
+
+        OUTPUT PARAMETERS:
+            RA      -   array of real numbers which stores model,
+                        array[0..RLen-1]
+            RLen    -   RA lenght
+
+          -- ALGLIB --
+             Copyright 15.03.2009 by Bochkanov Sergey
+        *************************************************************************/
+        public static void lrserialize(ref linearmodel lm,
+            ref double[] ra,
+            ref int rlen)
+        {
+            int i_ = 0;
+            int i1_ = 0;
+
+            rlen = (int)Math.Round(lm.w[0])+1;
+            ra = new double[rlen-1+1];
+            ra[0] = lrvnum;
+            i1_ = (0) - (1);
+            for(i_=1; i_<=rlen-1;i_++)
+            {
+                ra[i_] = lm.w[i_+i1_];
+            }
+        }
+
+
+        /*************************************************************************
+        Unserialization of DecisionForest strucure
+
+        INPUT PARAMETERS:
+            RA      -   real array which stores decision forest
+
+        OUTPUT PARAMETERS:
+            LM      -   unserialized structure
+
+          -- ALGLIB --
+             Copyright 15.03.2009 by Bochkanov Sergey
+        *************************************************************************/
+        public static void lrunserialize(ref double[] ra,
+            ref linearmodel lm)
+        {
+            int i_ = 0;
+            int i1_ = 0;
+
+            System.Diagnostics.Debug.Assert((int)Math.Round(ra[0])==lrvnum, "LRUnserialize: incorrect array!");
+            lm.w = new double[(int)Math.Round(ra[1])-1+1];
+            i1_ = (1) - (0);
+            for(i_=0; i_<=(int)Math.Round(ra[1])-1;i_++)
+            {
+                lm.w[i_] = ra[i_+i1_];
+            }
+        }
+
+
+        /*************************************************************************
+        Obsolete subroutine, use LRBuildS
+
+          -- ALGLIB --
+             Copyright 26.04.2008 by Bochkanov Sergey
+
+        References:
+        1. Numerical Recipes in C, "15.2 Fitting Data to a Straight Line"
+        *************************************************************************/
+        public static void lrlines(ref double[,] xy,
+            ref double[] s,
+            int n,
+            ref int info,
+            ref double a,
+            ref double b,
+            ref double vara,
+            ref double varb,
+            ref double covab,
+            ref double corrab,
+            ref double p)
+        {
+            int i = 0;
+            double ss = 0;
+            double sx = 0;
+            double sxx = 0;
+            double sy = 0;
+            double stt = 0;
+            double e1 = 0;
+            double e2 = 0;
+            double t = 0;
+            double chi2 = 0;
+
+            if( n<2 )
+            {
+                info = -1;
+                return;
+            }
+            for(i=0; i<=n-1; i++)
+            {
+                if( s[i]<=0 )
+                {
+                    info = -2;
+                    return;
+                }
+            }
+            info = 1;
+            
+            //
+            // Calculate S, SX, SY, SXX
+            //
+            ss = 0;
+            sx = 0;
+            sy = 0;
+            sxx = 0;
+            for(i=0; i<=n-1; i++)
+            {
+                t = AP.Math.Sqr(s[i]);
+                ss = ss+1/t;
+                sx = sx+xy[i,0]/t;
+                sy = sy+xy[i,1]/t;
+                sxx = sxx+AP.Math.Sqr(xy[i,0])/t;
+            }
+            
+            //
+            // Test for condition number
+            //
+            t = Math.Sqrt(4*AP.Math.Sqr(sx)+AP.Math.Sqr(ss-sxx));
+            e1 = 0.5*(ss+sxx+t);
+            e2 = 0.5*(ss+sxx-t);
+            if( Math.Min(e1, e2)<=1000*AP.Math.MachineEpsilon*Math.Max(e1, e2) )
+            {
+                info = -3;
+                return;
+            }
+            
+            //
+            // Calculate A, B
+            //
+            a = 0;
+            b = 0;
+            stt = 0;
+            for(i=0; i<=n-1; i++)
+            {
+                t = (xy[i,0]-sx/ss)/s[i];
+                b = b+t*xy[i,1]/s[i];
+                stt = stt+AP.Math.Sqr(t);
+            }
+            b = b/stt;
+            a = (sy-sx*b)/ss;
+            
+            //
+            // Calculate goodness-of-fit
+            //
+            if( n>2 )
+            {
+                chi2 = 0;
+                for(i=0; i<=n-1; i++)
+                {
+                    chi2 = chi2+AP.Math.Sqr((xy[i,1]-a-b*xy[i,0])/s[i]);
+                }
+                p = igammaf.incompletegammac(((double)(n-2))/(double)(2), chi2/2);
+            }
+            else
+            {
+                p = 1;
+            }
+            
+            //
+            // Calculate other parameters
+            //
+            vara = (1+AP.Math.Sqr(sx)/(ss*stt))/ss;
+            varb = 1/stt;
+            covab = -(sx/(ss*stt));
+            corrab = covab/Math.Sqrt(vara*varb);
+        }
+
+
+        /*************************************************************************
+        Obsolete subroutine, use LRBuild
+
+          -- ALGLIB --
+             Copyright 02.08.2008 by Bochkanov Sergey
+        *************************************************************************/
+        public static void lrline(ref double[,] xy,
+            int n,
+            ref int info,
+            ref double a,
+            ref double b)
+        {
+            double[] s = new double[0];
+            int i = 0;
+            double vara = 0;
+            double varb = 0;
+            double covab = 0;
+            double corrab = 0;
+            double p = 0;
+
+            if( n<2 )
+            {
+                info = -1;
+                return;
+            }
+            s = new double[n-1+1];
+            for(i=0; i<=n-1; i++)
+            {
+                s[i] = 1;
+            }
+            lrlines(ref xy, ref s, n, ref info, ref a, ref b, ref vara, ref varb, ref covab, ref corrab, ref p);
+        }
+
+
+        /*************************************************************************
+        Internal linear regression subroutine
+        *************************************************************************/
+        private static void lrinternal(ref double[,] xy,
+            ref double[] s,
+            int npoints,
+            int nvars,
+            ref int info,
+            ref linearmodel lm,
+            ref lrreport ar)
+        {
+            double[,] a = new double[0,0];
+            double[,] u = new double[0,0];
+            double[,] vt = new double[0,0];
+            double[,] vm = new double[0,0];
+            double[,] xym = new double[0,0];
+            double[] b = new double[0];
+            double[] sv = new double[0];
+            double[] t = new double[0];
+            double[] svi = new double[0];
+            double[] work = new double[0];
+            int i = 0;
+            int j = 0;
+            int k = 0;
+            int ncv = 0;
+            int na = 0;
+            int nacv = 0;
+            double r = 0;
+            double p = 0;
+            double epstol = 0;
+            lrreport ar2 = new lrreport();
+            int offs = 0;
+            linearmodel tlm = new linearmodel();
+            int i_ = 0;
+            int i1_ = 0;
+
+            epstol = 1000;
+            
+            //
+            // Check for errors in data
+            //
+            if( npoints<nvars | nvars<1 )
+            {
+                info = -1;
+                return;
+            }
+            for(i=0; i<=npoints-1; i++)
+            {
+                if( s[i]<=0 )
+                {
+                    info = -2;
+                    return;
+                }
+            }
+            info = 1;
+            
+            //
+            // Create design matrix
+            //
+            a = new double[npoints-1+1, nvars-1+1];
+            b = new double[npoints-1+1];
+            for(i=0; i<=npoints-1; i++)
+            {
+                r = 1/s[i];
+                for(i_=0; i_<=nvars-1;i_++)
+                {
+                    a[i,i_] = r*xy[i,i_];
+                }
+                b[i] = xy[i,nvars]/s[i];
+            }
+            
+            //
+            // Allocate W:
+            // W[0]     array size
+            // W[1]     version number, 0
+            // W[2]     NVars (minus 1, to be compatible with external representation)
+            // W[3]     coefficients offset
+            //
+            lm.w = new double[4+nvars-1+1];
+            offs = 4;
+            lm.w[0] = 4+nvars;
+            lm.w[1] = lrvnum;
+            lm.w[2] = nvars-1;
+            lm.w[3] = offs;
+            
+            //
+            // Solve problem using SVD:
+            //
+            // 0. check for degeneracy (different types)
+            // 1. A = U*diag(sv)*V'
+            // 2. T = b'*U
+            // 3. w = SUM((T[i]/sv[i])*V[..,i])
+            // 4. cov(wi,wj) = SUM(Vji*Vjk/sv[i]^2,K=1..M)
+            //
+            // see $15.4 of "Numerical Recipes in C" for more information
+            //
+            t = new double[nvars-1+1];
+            svi = new double[nvars-1+1];
+            ar.c = new double[nvars-1+1, nvars-1+1];
+            vm = new double[nvars-1+1, nvars-1+1];
+            if( !svd.rmatrixsvd(a, npoints, nvars, 1, 1, 2, ref sv, ref u, ref vt) )
+            {
+                info = -4;
+                return;
+            }
+            if( sv[0]<=0 )
+            {
+                
+                //
+                // Degenerate case: zero design matrix.
+                //
+                for(i=offs; i<=offs+nvars-1; i++)
+                {
+                    lm.w[i] = 0;
+                }
+                ar.rmserror = lrrmserror(ref lm, ref xy, npoints);
+                ar.avgerror = lravgerror(ref lm, ref xy, npoints);
+                ar.avgrelerror = lravgrelerror(ref lm, ref xy, npoints);
+                ar.cvrmserror = ar.rmserror;
+                ar.cvavgerror = ar.avgerror;
+                ar.cvavgrelerror = ar.avgrelerror;
+                ar.ncvdefects = 0;
+                ar.cvdefects = new int[nvars-1+1];
+                ar.c = new double[nvars-1+1, nvars-1+1];
+                for(i=0; i<=nvars-1; i++)
+                {
+                    for(j=0; j<=nvars-1; j++)
+                    {
+                        ar.c[i,j] = 0;
+                    }
+                }
+                return;
+            }
+            if( sv[nvars-1]<=epstol*AP.Math.MachineEpsilon*sv[0] )
+            {
+                
+                //
+                // Degenerate case, non-zero design matrix.
+                //
+                // We can leave it and solve task in SVD least squares fashion.
+                // Solution and covariance matrix will be obtained correctly,
+                // but CV error estimates - will not. It is better to reduce
+                // it to non-degenerate task and to obtain correct CV estimates.
+                //
+                for(k=nvars; k>=1; k--)
+                {
+                    if( sv[k-1]>epstol*AP.Math.MachineEpsilon*sv[0] )
+                    {
+                        
+                        //
+                        // Reduce
+                        //
+                        xym = new double[npoints-1+1, k+1];
+                        for(i=0; i<=npoints-1; i++)
+                        {
+                            for(j=0; j<=k-1; j++)
+                            {
+                                r = 0.0;
+                                for(i_=0; i_<=nvars-1;i_++)
+                                {
+                                    r += xy[i,i_]*vt[j,i_];
+                                }
+                                xym[i,j] = r;
+                            }
+                            xym[i,k] = xy[i,nvars];
+                        }
+                        
+                        //
+                        // Solve
+                        //
+                        lrinternal(ref xym, ref s, npoints, k, ref info, ref tlm, ref ar2);
+                        if( info!=1 )
+                        {
+                            return;
+                        }
+                        
+                        //
+                        // Convert back to un-reduced format
+                        //
+                        for(j=0; j<=nvars-1; j++)
+                        {
+                            lm.w[offs+j] = 0;
+                        }
+                        for(j=0; j<=k-1; j++)
+                        {
+                            r = tlm.w[offs+j];
+                            i1_ = (0) - (offs);
+                            for(i_=offs; i_<=offs+nvars-1;i_++)
+                            {
+                                lm.w[i_] = lm.w[i_] + r*vt[j,i_+i1_];
+                            }
+                        }
+                        ar.rmserror = ar2.rmserror;
+                        ar.avgerror = ar2.avgerror;
+                        ar.avgrelerror = ar2.avgrelerror;
+                        ar.cvrmserror = ar2.cvrmserror;
+                        ar.cvavgerror = ar2.cvavgerror;
+                        ar.cvavgrelerror = ar2.cvavgrelerror;
+                        ar.ncvdefects = ar2.ncvdefects;
+                        ar.cvdefects = new int[nvars-1+1];
+                        for(j=0; j<=ar.ncvdefects-1; j++)
+                        {
+                            ar.cvdefects[j] = ar2.cvdefects[j];
+                        }
+                        ar.c = new double[nvars-1+1, nvars-1+1];
+                        work = new double[nvars+1];
+                        blas.matrixmatrixmultiply(ref ar2.c, 0, k-1, 0, k-1, false, ref vt, 0, k-1, 0, nvars-1, false, 1.0, ref vm, 0, k-1, 0, nvars-1, 0.0, ref work);
+                        blas.matrixmatrixmultiply(ref vt, 0, k-1, 0, nvars-1, true, ref vm, 0, k-1, 0, nvars-1, false, 1.0, ref ar.c, 0, nvars-1, 0, nvars-1, 0.0, ref work);
+                        return;
+                    }
+                }
+                info = -255;
+                return;
+            }
+            for(i=0; i<=nvars-1; i++)
+            {
+                if( sv[i]>epstol*AP.Math.MachineEpsilon*sv[0] )
+                {
+                    svi[i] = 1/sv[i];
+                }
+                else
+                {
+                    svi[i] = 0;
+                }
+            }
+            for(i=0; i<=nvars-1; i++)
+            {
+                t[i] = 0;
+            }
+            for(i=0; i<=npoints-1; i++)
+            {
+                r = b[i];
+                for(i_=0; i_<=nvars-1;i_++)
+                {
+                    t[i_] = t[i_] + r*u[i,i_];
+                }
+            }
+            for(i=0; i<=nvars-1; i++)
+            {
+                lm.w[offs+i] = 0;
+            }
+            for(i=0; i<=nvars-1; i++)
+            {
+                r = t[i]*svi[i];
+                i1_ = (0) - (offs);
+                for(i_=offs; i_<=offs+nvars-1;i_++)
+                {
+                    lm.w[i_] = lm.w[i_] + r*vt[i,i_+i1_];
+                }
+            }
+            for(j=0; j<=nvars-1; j++)
+            {
+                r = svi[j];
+                for(i_=0; i_<=nvars-1;i_++)
+                {
+                    vm[i_,j] = r*vt[j,i_];
+                }
+            }
+            for(i=0; i<=nvars-1; i++)
+            {
+                for(j=i; j<=nvars-1; j++)
+                {
+                    r = 0.0;
+                    for(i_=0; i_<=nvars-1;i_++)
+                    {
+                        r += vm[i,i_]*vm[j,i_];
+                    }
+                    ar.c[i,j] = r;
+                    ar.c[j,i] = r;
+                }
+            }
+            
+            //
+            // Leave-1-out cross-validation error.
+            //
+            // NOTATIONS:
+            // A            design matrix
+            // A*x = b      original linear least squares task
+            // U*S*V'       SVD of A
+            // ai           i-th row of the A
+            // bi           i-th element of the b
+            // xf           solution of the original LLS task
+            //
+            // Cross-validation error of i-th element from a sample is
+            // calculated using following formula:
+            //
+            //     ERRi = ai*xf - (ai*xf-bi*(ui*ui'))/(1-ui*ui')     (1)
+            //
+            // This formula can be derived from normal equations of the
+            // original task
+            //
+            //     (A'*A)x = A'*b                                    (2)
+            //
+            // by applying modification (zeroing out i-th row of A) to (2):
+            //
+            //     (A-ai)'*(A-ai) = (A-ai)'*b
+            //
+            // and using Sherman-Morrison formula for updating matrix inverse
+            //
+            // NOTE 1: b is not zeroed out since it is much simpler and
+            // does not influence final result.
+            //
+            // NOTE 2: some design matrices A have such ui that 1-ui*ui'=0.
+            // Formula (1) can't be applied for such cases and they are skipped
+            // from CV calculation (which distorts resulting CV estimate).
+            // But from the properties of U we can conclude that there can
+            // be no more than NVars such vectors. Usually
+            // NVars << NPoints, so in a normal case it only slightly
+            // influences result.
+            //
+            ncv = 0;
+            na = 0;
+            nacv = 0;
+            ar.rmserror = 0;
+            ar.avgerror = 0;
+            ar.avgrelerror = 0;
+            ar.cvrmserror = 0;
+            ar.cvavgerror = 0;
+            ar.cvavgrelerror = 0;
+            ar.ncvdefects = 0;
+            ar.cvdefects = new int[nvars-1+1];
+            for(i=0; i<=npoints-1; i++)
+            {
+                
+                //
+                // Error on a training set
+                //
+                i1_ = (offs)-(0);
+                r = 0.0;
+                for(i_=0; i_<=nvars-1;i_++)
+                {
+                    r += xy[i,i_]*lm.w[i_+i1_];
+                }
+                ar.rmserror = ar.rmserror+AP.Math.Sqr(r-xy[i,nvars]);
+                ar.avgerror = ar.avgerror+Math.Abs(r-xy[i,nvars]);
+                if( xy[i,nvars]!=0 )
+                {
+                    ar.avgrelerror = ar.avgrelerror+Math.Abs((r-xy[i,nvars])/xy[i,nvars]);
+                    na = na+1;
+                }
+                
+                //
+                // Error using fast leave-one-out cross-validation
+                //
+                p = 0.0;
+                for(i_=0; i_<=nvars-1;i_++)
+                {
+                    p += u[i,i_]*u[i,i_];
+                }
+                if( p>1-epstol*AP.Math.MachineEpsilon )
+                {
+                    ar.cvdefects[ar.ncvdefects] = i;
+                    ar.ncvdefects = ar.ncvdefects+1;
+                    continue;
+                }
+                r = s[i]*(r/s[i]-b[i]*p)/(1-p);
+                ar.cvrmserror = ar.cvrmserror+AP.Math.Sqr(r-xy[i,nvars]);
+                ar.cvavgerror = ar.cvavgerror+Math.Abs(r-xy[i,nvars]);
+                if( xy[i,nvars]!=0 )
+                {
+                    ar.cvavgrelerror = ar.cvavgrelerror+Math.Abs((r-xy[i,nvars])/xy[i,nvars]);
+                    nacv = nacv+1;
+                }
+                ncv = ncv+1;
+            }
+            if( ncv==0 )
+            {
+                
+                //
+                // Something strange: ALL ui are degenerate.
+                // Unexpected...
+                //
+                info = -255;
+                return;
+            }
+            ar.rmserror = Math.Sqrt(ar.rmserror/npoints);
+            ar.avgerror = ar.avgerror/npoints;
+            if( na!=0 )
+            {
+                ar.avgrelerror = ar.avgrelerror/na;
+            }
+            ar.cvrmserror = Math.Sqrt(ar.cvrmserror/ncv);
+            ar.cvavgerror = ar.cvavgerror/ncv;
+            if( nacv!=0 )
+            {
+                ar.cvavgrelerror = ar.cvavgrelerror/nacv;
+            }
+        }
+    }
+}
Index: /trunk/sources/ALGLIB/normaldistr.cs
===================================================================
--- /trunk/sources/ALGLIB/normaldistr.cs	(revision 2445)
+++ /trunk/sources/ALGLIB/normaldistr.cs	(revision 2445)
@@ -0,0 +1,367 @@
+/*************************************************************************
+Cephes Math Library Release 2.8:  June, 2000
+Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
+
+Contributors:
+    * Sergey Bochkanov (ALGLIB project). Translation from C to
+      pseudocode.
+
+See subroutines comments for additional copyrights.
+
+>>> SOURCE LICENSE >>>
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation (www.fsf.org); either version 2 of the 
+License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+A copy of the GNU General Public License is available at
+http://www.fsf.org/licensing/licenses
+
+>>> END OF LICENSE >>>
+*************************************************************************/
+
+using System;
+
+namespace alglib
+{
+    public class normaldistr
+    {
+        /*************************************************************************
+        Error function
+
+        The integral is
+
+                                  x
+                                   -
+                        2         | |          2
+          erf(x)  =  --------     |    exp( - t  ) dt.
+                     sqrt(pi)   | |
+                                 -
+                                  0
+
+        For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
+        erf(x) = 1 - erfc(x).
+
+
+        ACCURACY:
+
+                             Relative error:
+        arithmetic   domain     # trials      peak         rms
+           IEEE      0,1         30000       3.7e-16     1.0e-16
+
+        Cephes Math Library Release 2.8:  June, 2000
+        Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
+        *************************************************************************/
+        public static double erf(double x)
+        {
+            double result = 0;
+            double xsq = 0;
+            double s = 0;
+            double p = 0;
+            double q = 0;
+
+            s = Math.Sign(x);
+            x = Math.Abs(x);
+            if( x<0.5 )
+            {
+                xsq = x*x;
+                p = 0.007547728033418631287834;
+                p = 0.288805137207594084924010+xsq*p;
+                p = 14.3383842191748205576712+xsq*p;
+                p = 38.0140318123903008244444+xsq*p;
+                p = 3017.82788536507577809226+xsq*p;
+                p = 7404.07142710151470082064+xsq*p;
+                p = 80437.3630960840172832162+xsq*p;
+                q = 0.0;
+                q = 1.00000000000000000000000+xsq*q;
+                q = 38.0190713951939403753468+xsq*q;
+                q = 658.070155459240506326937+xsq*q;
+                q = 6379.60017324428279487120+xsq*q;
+                q = 34216.5257924628539769006+xsq*q;
+                q = 80437.3630960840172826266+xsq*q;
+                result = s*1.1283791670955125738961589031*x*p/q;
+                return result;
+            }
+            if( x>=10 )
+            {
+                result = s;
+                return result;
+            }
+            result = s*(1-erfc(x));
+            return result;
+        }
+
+
+        /*************************************************************************
+        Complementary error function
+
+         1 - erf(x) =
+
+                                  inf.
+                                    -
+                         2         | |          2
+          erfc(x)  =  --------     |    exp( - t  ) dt
+                      sqrt(pi)   | |
+                                  -
+                                   x
+
+
+        For small x, erfc(x) = 1 - erf(x); otherwise rational
+        approximations are computed.
+
+
+        ACCURACY:
+
+                             Relative error:
+        arithmetic   domain     # trials      peak         rms
+           IEEE      0,26.6417   30000       5.7e-14     1.5e-14
+
+        Cephes Math Library Release 2.8:  June, 2000
+        Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
+        *************************************************************************/
+        public static double erfc(double x)
+        {
+            double result = 0;
+            double p = 0;
+            double q = 0;
+
+            if( x<0 )
+            {
+                result = 2-erfc(-x);
+                return result;
+            }
+            if( x<0.5 )
+            {
+                result = 1.0-erf(x);
+                return result;
+            }
+            if( x>=10 )
+            {
+                result = 0;
+                return result;
+            }
+            p = 0.0;
+            p = 0.5641877825507397413087057563+x*p;
+            p = 9.675807882987265400604202961+x*p;
+            p = 77.08161730368428609781633646+x*p;
+            p = 368.5196154710010637133875746+x*p;
+            p = 1143.262070703886173606073338+x*p;
+            p = 2320.439590251635247384768711+x*p;
+            p = 2898.0293292167655611275846+x*p;
+            p = 1826.3348842295112592168999+x*p;
+            q = 1.0;
+            q = 17.14980943627607849376131193+x*q;
+            q = 137.1255960500622202878443578+x*q;
+            q = 661.7361207107653469211984771+x*q;
+            q = 2094.384367789539593790281779+x*q;
+            q = 4429.612803883682726711528526+x*q;
+            q = 6089.5424232724435504633068+x*q;
+            q = 4958.82756472114071495438422+x*q;
+            q = 1826.3348842295112595576438+x*q;
+            result = Math.Exp(-AP.Math.Sqr(x))*p/q;
+            return result;
+        }
+
+
+        /*************************************************************************
+        Normal distribution function
+
+        Returns the area under the Gaussian probability density
+        function, integrated from minus infinity to x:
+
+                                   x
+                                    -
+                          1        | |          2
+           ndtr(x)  = ---------    |    exp( - t /2 ) dt
+                      sqrt(2pi)  | |
+                                  -
+                                 -inf.
+
+                    =  ( 1 + erf(z) ) / 2
+                    =  erfc(z) / 2
+
+        where z = x/sqrt(2). Computation is via the functions
+        erf and erfc.
+
+
+        ACCURACY:
+
+                             Relative error:
+        arithmetic   domain     # trials      peak         rms
+           IEEE     -13,0        30000       3.4e-14     6.7e-15
+
+        Cephes Math Library Release 2.8:  June, 2000
+        Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
+        *************************************************************************/
+        public static double normaldistribution(double x)
+        {
+            double result = 0;
+
+            result = 0.5*(erf(x/1.41421356237309504880)+1);
+            return result;
+        }
+
+
+        /*************************************************************************
+        Inverse of the error function
+
+        Cephes Math Library Release 2.8:  June, 2000
+        Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
+        *************************************************************************/
+        public static double inverf(double e)
+        {
+            double result = 0;
+
+            result = invnormaldistribution(0.5*(e+1))/Math.Sqrt(2);
+            return result;
+        }
+
+
+        /*************************************************************************
+        Inverse of Normal distribution function
+
+        Returns the argument, x, for which the area under the
+        Gaussian probability density function (integrated from
+        minus infinity to x) is equal to y.
+
+
+        For small arguments 0 < y < exp(-2), the program computes
+        z = sqrt( -2.0 * log(y) );  then the approximation is
+        x = z - log(z)/z  - (1/z) P(1/z) / Q(1/z).
+        There are two rational functions P/Q, one for 0 < y < exp(-32)
+        and the other for y up to exp(-2).  For larger arguments,
+        w = y - 0.5, and  x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
+
+        ACCURACY:
+
+                             Relative error:
+        arithmetic   domain        # trials      peak         rms
+           IEEE     0.125, 1        20000       7.2e-16     1.3e-16
+           IEEE     3e-308, 0.135   50000       4.6e-16     9.8e-17
+
+        Cephes Math Library Release 2.8:  June, 2000
+        Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
+        *************************************************************************/
+        public static double invnormaldistribution(double y0)
+        {
+            double result = 0;
+            double expm2 = 0;
+            double s2pi = 0;
+            double x = 0;
+            double y = 0;
+            double z = 0;
+            double y2 = 0;
+            double x0 = 0;
+            double x1 = 0;
+            int code = 0;
+            double p0 = 0;
+            double q0 = 0;
+            double p1 = 0;
+            double q1 = 0;
+            double p2 = 0;
+            double q2 = 0;
+
+            expm2 = 0.13533528323661269189;
+            s2pi = 2.50662827463100050242;
+            if( y0<=0 )
+            {
+                result = -AP.Math.MaxRealNumber;
+                return result;
+            }
+            if( y0>=1 )
+            {
+                result = AP.Math.MaxRealNumber;
+                return result;
+            }
+            code = 1;
+            y = y0;
+            if( y>1.0-expm2 )
+            {
+                y = 1.0-y;
+                code = 0;
+            }
+            if( y>expm2 )
+            {
+                y = y-0.5;
+                y2 = y*y;
+                p0 = -59.9633501014107895267;
+                p0 = 98.0010754185999661536+y2*p0;
+                p0 = -56.6762857469070293439+y2*p0;
+                p0 = 13.9312609387279679503+y2*p0;
+                p0 = -1.23916583867381258016+y2*p0;
+                q0 = 1;
+                q0 = 1.95448858338141759834+y2*q0;
+                q0 = 4.67627912898881538453+y2*q0;
+                q0 = 86.3602421390890590575+y2*q0;
+                q0 = -225.462687854119370527+y2*q0;
+                q0 = 200.260212380060660359+y2*q0;
+                q0 = -82.0372256168333339912+y2*q0;
+                q0 = 15.9056225126211695515+y2*q0;
+                q0 = -1.18331621121330003142+y2*q0;
+                x = y+y*y2*p0/q0;
+                x = x*s2pi;
+                result = x;
+                return result;
+            }
+            x = Math.Sqrt(-(2.0*Math.Log(y)));
+            x0 = x-Math.Log(x)/x;
+            z = 1.0/x;
+            if( x<8.0 )
+            {
+                p1 = 4.05544892305962419923;
+                p1 = 31.5251094599893866154+z*p1;
+                p1 = 57.1628192246421288162+z*p1;
+                p1 = 44.0805073893200834700+z*p1;
+                p1 = 14.6849561928858024014+z*p1;
+                p1 = 2.18663306850790267539+z*p1;
+                p1 = -(1.40256079171354495875*0.1)+z*p1;
+                p1 = -(3.50424626827848203418*0.01)+z*p1;
+                p1 = -(8.57456785154685413611*0.0001)+z*p1;
+                q1 = 1;
+                q1 = 15.7799883256466749731+z*q1;
+                q1 = 45.3907635128879210584+z*q1;
+                q1 = 41.3172038254672030440+z*q1;
+                q1 = 15.0425385692907503408+z*q1;
+                q1 = 2.50464946208309415979+z*q1;
+                q1 = -(1.42182922854787788574*0.1)+z*q1;
+                q1 = -(3.80806407691578277194*0.01)+z*q1;
+                q1 = -(9.33259480895457427372*0.0001)+z*q1;
+                x1 = z*p1/q1;
+            }
+            else
+            {
+                p2 = 3.23774891776946035970;
+                p2 = 6.91522889068984211695+z*p2;
+                p2 = 3.93881025292474443415+z*p2;
+                p2 = 1.33303460815807542389+z*p2;
+                p2 = 2.01485389549179081538*0.1+z*p2;
+                p2 = 1.23716634817820021358*0.01+z*p2;
+                p2 = 3.01581553508235416007*0.0001+z*p2;
+                p2 = 2.65806974686737550832*0.000001+z*p2;
+                p2 = 6.23974539184983293730*0.000000001+z*p2;
+                q2 = 1;
+                q2 = 6.02427039364742014255+z*q2;
+                q2 = 3.67983563856160859403+z*q2;
+                q2 = 1.37702099489081330271+z*q2;
+                q2 = 2.16236993594496635890*0.1+z*q2;
+                q2 = 1.34204006088543189037*0.01+z*q2;
+                q2 = 3.28014464682127739104*0.0001+z*q2;
+                q2 = 2.89247864745380683936*0.000001+z*q2;
+                q2 = 6.79019408009981274425*0.000000001+z*q2;
+                x1 = z*p2/q2;
+            }
+            x = x0-x1;
+            if( code!=0 )
+            {
+                x = -x;
+            }
+            result = x;
+            return result;
+        }
+    }
+}
Index: /trunk/sources/HeuristicLab.LinearRegression/3.2/LinearRegressionOperator.cs
===================================================================
--- /trunk/sources/HeuristicLab.LinearRegression/3.2/LinearRegressionOperator.cs	(revision 2444)
+++ /trunk/sources/HeuristicLab.LinearRegression/3.2/LinearRegressionOperator.cs	(revision 2445)
@@ -101,10 +101,24 @@
 
     private double[] CalculateCoefficients(double[,] inputMatrix, double[] targetVector) {
-      double[] weights = new double[targetVector.Length];
-      double[] coefficients = new double[inputMatrix.GetLength(1)];
-      for (int i = 0; i < weights.Length; i++) weights[i] = 1.0;
-      // call external ALGLIB solver
-      alglib.leastsquares.buildgeneralleastsquares(ref targetVector, ref weights, ref inputMatrix, inputMatrix.GetLength(0), inputMatrix.GetLength(1), ref coefficients);
+      int retVal = 0;
+      alglib.linreg.linearmodel lm = new alglib.linreg.linearmodel();
+      alglib.linreg.lrreport ar = new alglib.linreg.lrreport();
+      int n = targetVector.Length;
+      int p = inputMatrix.GetLength(1);
+      double[,] dataset = new double[n, p];
+      for (int row = 0; row < n; row++) {
+        for (int column = 0; column < p-1; column++) {
+          dataset[row, column] = inputMatrix[row, column];
+        }
+        dataset[row, p-1] = targetVector[row];
+      }
+      alglib.linreg.lrbuild(ref dataset, n, p-1, ref retVal, ref lm, ref ar);
+      if (retVal != 1) throw new ArgumentException("Error in calculation of linear regression model");
+      Console.Out.WriteLine("ALGLIB Linear Regression: Estimated generalization RMS = {0}", ar.cvrmserror);
 
+      double[] coefficients = new double[p];
+      for (int i = 0; i < p; i++) {
+        coefficients[i] = lm.w[i+4];
+      }
       return coefficients;
     }
