/**********************************************************/
/*   Program Regresion_Polinomial_Por_Minimos_Cuadrados   */
/**********************************************************/

#include <stdio.h>
#include <conio.h>
#include <math.h>

void Lee_Datos(void);
void Calcula_Matriz_Del_Sistema(void);
void Imprime_Sistema(void);
void Imprime_Linea_Matriz(void);
void Eliminacion_Gaussiana(void);
void Error_Por_Minimos_Cuadrados(void);

  float  A[12][12];
  float  Fx[12];
  float  X[12];
  float  Xr[12];
  int    N,M;
  int    k,j,i;
  int    Nx;
  float  Pot;
  float  Error, Error_med;
  float  Sum_Error2, Sum_ErrorM, SumY;
  float  Fx_Cal, Ymed;
  float  Resul;
  float  Piv;
  float  Sum;

/**********   PROGRAMA PRINCIPAL   **********/

int main(void)
{
	 Lee_Datos();
	 Calcula_Matriz_Del_Sistema();
	 Imprime_Sistema();
	 Eliminacion_Gaussiana();
	 Error_Por_Minimos_Cuadrados();
}

void Lee_Datos(void)

  {
	 clrscr();
	 printf("\n  Numero de Datos =");
	 scanf("%d",&N);
	 printf("\n  Grado del Polinomio =");
	 scanf("%d",&M);

	 printf("\n       Valores Tabulados ");
	 printf("\n--------------------------------");
	 printf("\n i       X           Fx");
	 printf("\n--------------------------------");

	 X[0] = 0;      Fx[0]= 10;
	 X[1] = 2;     Fx[1]= 12;
	 X[2] = 4;     Fx[2]= 18;
	 X[3] = 4;     Fx[3]= 22;
	 X[4] = 8;     Fx[4]= 20;
	 X[5] = 12;     Fx[5]= 30;
	 X[6] = 16;     Fx[6]= 26;
	 X[7] = 20;      Fx[7]= 30;
	 X[8] = 24;     Fx[8]= 26;
	 X[9] = 28;     Fx[9]= 28;
	 X[10] = 30;    Fx[10]= 22;
	 X[11] = 34;    Fx[11]= 20;

	 SumY = 0;
	 for (k=0; k<N;k++)
		{
		  printf("\n%d %10.4f  |  %10.4f",k+1,X[k],Fx[k]);
			SumY += Fx[k];
		}
	 Ymed = SumY/N;
	 printf("\n--------------------------------");
  }

void Calcula_Matriz_Del_Sistema(void)

	/*  Calculo de las sumatorias para formar el sistema */

{
	for (k=0; k<M+1; k++ )
    {
		for (j=0;j<M+2;j++)
        {
			 if ( k+j==0 )
				 A[k][j] = N;
			 else
				 {
					Pot =k+j;
					A[k][j]=0;
					for (i=0; i<N; i++)
                  {
						  if (j == M+1)
							  A[k][j] += pow(X[i],k)*Fx[i];
						  else
							  A[k][j] += pow(X[i],Pot);
                  }
              }
        }
	 }
}

void Imprime_Sistema(void)

{
	 printf("\n Sistema de Ecuaciones Lineales:");
	 printf("\n--------------------------------\n");

	 for (k=0; k< M+1; k++)
      {
		  for (j=0; j<=M+1; j++)
          {
				if (j == M+1)
					printf("  =  %10.4f",A[k][j]);
				else
					printf("%10.4f a%d",A[k][j],j);
          }
		  printf("\n");
      }
	 printf("\n");
getche();
}


void Imprime_Linea_Matriz(void)

  {
	  for (j=0; j<=M+1;j++)
       {
			if (j == M+1)
				printf("  =  %10.4f",A[k][j]);
			else
				printf("%10.4f a%d",A[k][j],j);
		 }
	}

void Eliminacion_Gaussiana(void)

{
  printf("\n-------------ELIMINACION HACIA ADELANTE --------------");
  printf("\n            Ecuacion Triangula Superior:\n");
  for (k=0;k< M+1;k++)
   {
	  for (i=k+1;i<M+1;i++)
       {
			Piv=A[i][k]/A[k][k];
			for (j=k+1;j<=M+1;j++)
           {
				 A[i][j] = A[i][j] - Piv*A[k][j];
			  }
		 }
	  for (i=k+1;i<M+1;i++)
       {
			A[i][k] =0;
		 }
	  Imprime_Linea_Matriz();
	  printf("\n");
	}
  printf("\n------------------------------------------------------");
  printf( "\n------------------ELIMINACION HACIA ATRAS------------");
  printf("\n------------------------------------------------------");
  for (i = 0; i < M+1; i++)
    {
      Sum = 0;
		k = M - i;
		for (j=(M-i); j < M+1; j++)
        {
			Sum += A[k][j] * Xr[j];
        }
		Xr[k] = ( A[k][M+1] - Sum ) / A[k][k];
		printf("\nLa a(%d) es igua a =>%10.4f",k,Xr[k]);
	 }
	 getche();
}

void Error_Por_Minimos_Cuadrados(void)

/*   Calculo de Errores Al cuadrado   */

{
	Sum_Error2=0;
	Sum_ErrorM=0;
	printf("\n Valores Estadisticos del Modelo");
	printf("\n-----------------------------------------------");
	printf("\n i      Xi     Fxi       (Yi-Ymed)^2    (Sr)^2");
	printf("\n------------------------------------------------");
	for (i=0; i< N; i++)
     {
		 Error = 0;
		 Fx_Cal = 0;
		 Error_med = pow((Fx[i]-Ymed),2);
		 for (j=0; j<M+1; j++)
          {
				if (j==0)
					Fx_Cal +=  Xr[j];
				else
					Fx_Cal +=  pow(X[i],j)*Xr[j];
			 }
       Error = (Fx[i] - Fx_Cal)*(Fx[i] - Fx_Cal);
		 printf("\n %2d%10.2f%10.2f%12.4f%12.4f",i+1,X[i],Fx[i],Error_med,Error);
		 Sum_Error2 +=  Error;
		 Sum_ErrorM += Error_med;
	  }
	printf("\n------------------------------------------------");
	printf("\nSumatoria de los Errores al Cuadrados..(Sr)  => %8.6f",Sum_Error2);
	printf("\nSumatoria de los Errores a la media....(St)  => %8.6f",Sum_ErrorM);
	printf("\nCalculo del Desviacion estandar........(Sy)  => %8.6f",sqrt(Sum_ErrorM/(N-1)));
	printf("\nCalculo del Error estandar.............(Sy/x)=> %8.6f",sqrt(Sum_Error2/(N-(M+1))));
	printf("\nCoeficiente de determinacion...........(r^2) => %8.6f",(Sum_ErrorM-Sum_Error2)/(Sum_ErrorM));
	printf("\nCoeficiente de correlacion.............(r)   => %8.6f",sqrt((Sum_ErrorM-Sum_Error2)/(Sum_ErrorM)));

	getche();
}