RÉGRESSION LINÉAIRE MULTIPLE: PLAN DES MOINDRES CARRÉS

// regressionmultiple.cpp : Defines the entry point for the console application.
//

// regressionmultiple.cpp : Defines the entry point for the console application.
//

#include <stdio.h> 
#include <stdlib.h> 
#include <time.h> 

/* Nombre de points 3D dans le nuage de points */ 
#define NUM_POINT 20 

/* Valeur maximale d'un nombre aléatoire généré pour X ,Y et Z (100 = un nombre aléatoire entre 0 et 99) */ 
static const unsigned int RX = 100; 
static const unsigned int RY = 100; 
static const unsigned int RZ = 100; 

typedef struct 
{ 
  int x; 
  int y; 
  int z; 
} Point3D; 

typedef struct 
{ 
  float a; 
  float b; 
  float c; 
} Res; 

void generateRandomPoints(Point3D vect[]) 
{ 
  unsigned int i; 
  for (i = 0; i < NUM_POINT; ++i) 
  { 
    vect[i].x = rand() % RX; 
    vect[i].y = rand() % RY; 
    vect[i].z = rand() % RZ; 
  } 
} 

Res multipleLinearRegression(Point3D vect[]) 
{ 
  float sumxi2 = 0; 
  float sumxiyi = 0; 
  float sumxi = 0; 
  float sumyi2 = 0; 
  float sumyi = 0; 
  float sumxizi = 0; 
  float sumyizi = 0; 
  float sumzi = 0; 
  Point3D pt; 
  unsigned int i; 
  Res res; 

  /* Regression multiple selon la methode des moindres carrés */ 
  for (i = 0; i < NUM_POINT; ++i) 
  { 
    pt = vect[i]; 
    sumxi2 += pt.x * pt.x; 
    sumxiyi += pt.x * pt.y; 
    sumxi += pt.x; 
    sumyi2 += pt.y * pt.y; 
    sumyi += pt.y; 
    sumxizi += pt.x * pt.z; 
    sumyizi += pt.y * pt.z; 
    sumzi += pt.z; 
  } 

  /* Equation du plan : ax + by + c = z */ 
  res.c = ((sumxi * sumxizi - sumxi2 * sumzi) * (-sumxiyi * sumyi + sumxi * sumyi2) - 
    ((-sumxi2 * sumyi + sumxi * sumxiyi) * (sumxi * sumyizi - sumxiyi * sumzi))) / 
    (-((-sumxi2 * sumyi + sumxi * sumxiyi) * (-sumxiyi * NUM_POINT + sumxi * sumyi)) + 
     (-sumxiyi * sumyi + sumxi * sumyi2) * (-sumxi2 * NUM_POINT + sumxi * sumxi)); 
  res.b = ((sumxi * sumyizi - sumxiyi * sumzi) - (-sumxiyi * NUM_POINT + sumxi * sumyi) * res.c) / 
    (-sumxiyi * sumyi + sumxi * sumyi2); 
  res.a = (sumzi - sumyi * res.b - NUM_POINT * res.c) / sumxi; 

  return res; 
} 

void writeCloud(FILE* file, const Point3D vect[]) 
{ 
  unsigned int i; 

  for (i = 0; i < NUM_POINT; ++i) 
    fprintf(file, "%i\t%i\t%i\n", vect[i].x, vect[i].y, vect[i].z); 
} 

void writePlan(FILE* file, const Res res) 
{ 
  float x, y; 

  for (x = 0; x < RX; ++x) 
    for (y = 0; y < RY; ++y) 
      fprintf(file, "%f\t%f\t%f\n", x, x, res.a * x + res.b * y + res.c); 
} 

void writeCombination(FILE* file, const Point3D vect[], const Res res) 
{ 
  writeCloud(file, vect); 
  writePlan(file, res); 
} 

void writeToFile(const char* fileplan, const char* filecloud, const char* filecombined, const Point3D vect[], const Res res) 
{ 
  /* Ecriture dans un fichier des points qui peuvent etre affiché dans un editeur graphique tel que gnuplot, excel, ect. */ 
  /* Le dossier tmp doit exister, dans le cas contraire, créer le dossier tmp dans à la racine c: */ 
  /* Les points sont formatés par ligne avec comme valeur X, la premiere colonne, Y pour la deuxieme et Z pour la troisieme. */ 
  FILE* file = NULL; 

  file = fopen(fileplan, "w");/* "c:\\tmp\\nuagedepoint.txt" */ 
  writeCloud(file, vect); 
  fclose(file); 

  file = fopen(filecloud, "w"); /* "c:\\tmp\\plan.txt" */ 
  writePlan(file, res); 
  fclose(file); 

  file = fopen(filecombined, "w"); /* "c:\\tmp\\combinaisonNuagePlan.txt" */ 
  writeCombination(file, vect, res); 
  fclose(file); 
} 

int main(void) 
{ 
  Point3D vect[NUM_POINT]; 
  Res res; 
  srand(time(NULL)); 

  res = multipleLinearRegression(vect); 
  printf("Le plan minimisant les erreurs sur le nuages de points est : %f x + %f y + %f = z\n", res.a, res.b, res.c); 

  writeToFile("/tmp/plan", "/tmp/cloud", "/tmp/combined", vect, res); 

  /* getchar(); */ 
  /* system("pause"); */ 

  return 0; 
}