using System;

using Bluebit.MatrixLibrary;



namespace Regression

{

    // This example demonstrates how to generate a polynomial regression

    // for a given data set of x,y pairs.

    //

    // A 3rd order polynomial is used.

    // y = a + b*x + c*x^2 + d*x^3

    class Test

    {

        [STAThread]

	static void Main(string[] args)

	{

            // Array X contains the values of of x

            double[] xArray = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16};

            // array Y contains the observed values of y

            double[] yArray = {1,3,2,4,4,6,3,5,7,8, 7,10, 8, 9,10,12,15};



            // order of the polynomial

            int n = 3;



            // number or (X,Y) pairs

            int count = xArray.Length;



            // Create the X matrix containg the powers of X values

            Matrix X = new Matrix(count,n+1); 

            for(int i=0; i< count ;i++)

                for(int j=0; j<=n; j++)

                    X[i,j] = Math.Pow(xArray[i], (double) j);





            // Create the Y matrix containg the values of Y

            Matrix Y = new Matrix(yArray, count, 1, StorageOrder.ColumnMajor) ;



            // Use QR decomposition of matrix X to perform a least squares fit

            QR qr = new QR(X);



            // Matrix S contains the solution

            Matrix S = qr.Solve(Y);



            // Matrix Yest contains the estimated values of Y

            Matrix Yest = X * S;



            // Matrix dY will contains the differences Y - Yest;

            Matrix dY = Y - Yest;



            // The sumSq is the sum of squares 

            double sumSq = dY.ColsDotProduct(0,0);  



            // print the results

            Console.WriteLine("Polynomial regression y = a + b*x + c*x^2 + d*x^3 ");

            Console.WriteLine("Calculated coefficients of the polynomial:");

            

            Console.WriteLine(" a = {0,9:N6}",S[0,0]);

            Console.WriteLine(" b = {0,9:N6}",S[1,0]);

            Console.WriteLine(" c = {0,9:N6}",S[2,0]);     

            Console.WriteLine(" d = {0,9:N6}",S[3,0]);

            Console.WriteLine();

            Console.Write("y = ");

            for(int i = 0; i <= n; i++)

                Console.Write(" {0:+ 0.000;- 0.000}*x^{1}", S[i,0], i );

            Console.WriteLine("");

          

            Console.WriteLine("     x value y observed  y estimated   ");  

            for(int i=0; i<count;i++)

            Console.WriteLine(" {0,10:F3}  {1,10:F3}  {2,10:F3}", X[i,1], Y[i,0], Yest[i,0]);



            Console.Read();

        }

    }

}