Défi2 = les moindres carrés

' ******************************************************************
' Generateur de nombres aleatoires gaussiens
' ******************************************************************

DIM GaussSave  ' Sauvegarde un nombre aleatoire
DIM GaussNew%  ' Marque un nouveau calcul

FUNCTION RanGauss()
' ------------------------------------------------------------------
' Algorithme de Box-Muller :
' Calcule 2 nombre aleatoires de la loi normale reduite
' La fonction retourne un nombre et sauvegarde l'autre
' pour l'appel suivant
' ------------------------------------------------------------------

  DIM X, R, Theta

  IF GaussNew = 0 THEN
    X = RND() : R = SQR(-2 * LOG(X))
    X = RND() : Theta = 6.283185307179586 * X

    RanGauss = R * COS(Theta)  ' Retourne le premier nombre
    GaussSave = R * SIN(Theta) ' Sauvegarde le second
  ELSE
    RanGauss = GaussSave       ' Retourne le nombre sauvegardé
  END IF

  GaussNew = NOT GaussNew
END_FUNCTION

? RanGauss()

rem ============================================================================
rem                  Ajustement statistique                      
rem       Régression par la méthode des moindres carrés
rem                Papydall le 24 / 10 / 2019
rem ============================================================================
rem Le programme ci-dessous vous fournira les caractéristiques de l’ajustement
rem linéaire de type f(x) = ax + b d’un nuage de points Mi(xi,yi) en minimisant
rem la somme des carrées des distances entre yi et f(xi).
rem Le programme fournira aussi le point moyen, les variances de X et Y et leur
rem écarts-types, la covariance du couple (X,Y) et le coefficient de corrélation.
rem Le programme vous fournira également les valeurs des cinq sommes fondamentales
rem nécessaires à l’obtention de ces valeurs
rem ============================================================================

dim x(20), y(20) : ' coordonnées des points
dim s(10), w(10)
dim a(12,12)
dim np : ' nombre de points
dim p  : ' degré du lissage, degré du polynome
dim sx, sy, sxy, sx2, sy2
dim mx, my, mxy, vx, vy, cov, ect_x, ect_y
dim v, somme, a1, a2, b1, b2, nmi, corr
dim wg, hg, unite, xo,yo, npt
rem ============================================================================
width 0,800
picture 10 : full_space 10 : 2d_target_is 10 : print_target_is 10 : font_bold 10
alpha 20 : left 20,500 : top 20,100 : font_bold 20 : color 20,255,255,255
caption 0,"Ajustement statistique : Régression par la méthode des moindres carrés"
wg = 400 : hg = 400 : xo = 50 : yo = 420 : p = 1
rem ============================================================================
data 16 : ' Nombre de points, pour cet exemple 16 points
data 1,0.5, 2,1.8, 3,1.8, 4,2.2, 5,2.7, 6,3.3,  7,3.5, 8.2,4
data 9,4.7, 10,5,  11,6,  12,6,  13,7,  14,7.8, 15,9,  16,9.5
rem ============================================================================
Init_Points() : unite = wg / (np+1)
Grille()
Nuage()
Calcul()
Tracer()
end
rem ============================================================================
SUB Init_Points()
   dim_local i
   read np : npt = np+1
   for i = 1 to np : read x(i) : read y(i) : next i
END_SUB
rem ============================================================================
SUB Grille()
   dim_local i,x
   2d_line xo,yo,xo,5               : ' Axe des Y
   2d_line xo-10,yo-10,wg+100,hg+10 : ' Axe des X
   print_locate xo-10,wg+20 : print 0
   print_locate wg+100,hg+3 : print ">"
   print_locate xo-3,1      : print "^"
   2d_pen_dot : ' tracé avec séries de points
   for i = xo+unite to wg+80 step unite
       2d_line i,hg+10,i,10 : x = x+1
       print_locate i-5,hg+20 : print x
   next i
   x = 0
   for i = hg+10-unite to 10 step 0-unite
       2d_line xo,i,wg+80,i : x = x+1
       print_locate xo-20,i-5 : print x
   next i
   2d_pen_solid : ' tracé avec trait plein
END_SUB
rem ============================================================================
SUB Nuage()
   dim_local xp,yp,i    
   2d_fill_color 255,0,0
   for i = 1 to np
       xp = x(i)*unite : yp = y(i)*unite : ' Adaptation aux coordonnées de la grille
       2d_circle xo+xp,hg+10-yp,4        : ' Marquer le point
   next i
END_SUB
rem ============================================================================
SUB Calcul()
   dim_local i,j,k,t$
   for i = 1 to np
       sx  = sx + x(i)       : sy  = sy + y(i) : sxy = sxy + x(i)*y(i)
       sx2 = sx2 + x(i)*x(i) : sy2 = sy2 + y(i)*y(i)
   next i
   mx = sx/np : my = sy/np : mxy = sxy/np
   vx = sx2/np - mx*mx : vy = sy2/np-my*my
   cov = mxy-mx*my : ect_x = sqr(vx) : ect_y = sqr(vy)
   
   for k = 0 to 2*p
       s(k) = 0
       for i = 1 to np : s(k) = s(k) + power(x(i),k) : next i
   next k
   for k = 0 to p
       w(k) = 0
       for i = 1 to np : w(k) = w(k) + y(i)*power(x(i),k) : next i
   next k
   
   for k = 0 to p : a(k+1,p+2) = w(k) : next k
   for i = 1 to p+1
       for j = 1 to p+1 : a(i,j) = s(i+j-2) : next j
   next i
   a(1,1) = np  
   np = p+1           : ' dimension du système des équations
   for k = 1 to np-1  : ' résolution du système par la méthode du pivot
       v = a(k,k)
       if v = 0
          v = permuter(k)
          if v = 0 then singulier() : exit_sub
       end_if    
       for i = k+1 to np
           for j = k+1 to np+1
               a(i,j) = a(i,j) - a(i,k)*a(k,j)/v
           next j
       next i
   next k
   if a(np,np) = 0 then singulier() : exit_sub
   x(np) = a(np,np+1) / a(np,np)
   for i = np-1 to 1 step -1
       somme = 0
       for j = i+1 to np
           somme = somme + a(i,j) * x(j)
       next j
       x(i) = (a(i,np+1) - somme) / a(i,i)  
   next i
   for k = 1 to np : y(k) = x(np-k+1) : next k
   for k = 1 to np : x(k) = y(k) : next k
   a1 = cov/vx : a2 = cov/vy : b1 = my-a1*mx : b2 = mx-a2*my
   corr = cov/sqr(vx)/sqr(vy)  
   t$ = "La droite y = ax + b : a = " + str$(int(a1*1E3)/1E3)  
   t$ = t$ + "   b = " + str$(int(b1*1E3)/1E3) + chr$(13)  
   t$ = t$ + "Somme des xi : Sxi = " + str$(int(sx*1E3)/1E3) +chr$(13)
   t$ = t$ + "Somme des yi : Syi = " + str$(int(sy*1E3)/1E3) + chr$(13)
   t$ = t$ + "Somme des xi carrés : Sxi² = " + str$(int(sx2*1E3)/1E3) + chr$(13)
   t$ = t$ + "Somme des yi carrés : Syi² = " + str$(int(sy2*1E3)/1E3) + chr$(13)
   t$ = t$ + "Somme du produit des xi par yi : Sxiyi = " + str$(int(sxy*1E3)/1E3)+ chr$(13)    
   t$ = t$ + "Moyenne des X : Xmoyen = " + str$(int(mx*1E3)/1E3) +chr$(13)
   t$ = t$ + "Moyenne des Y : Ymoyen = " + str$(int(my*1E3)/1E3) + chr$(13)
   t$ = t$ + "Variance de X : V(x) = " + str$(int(vx*1E3)/1E3) +chr$(13)
   t$ = t$ + "Variance de Y : V(y) = " + str$(int(vy*1E3)/1E3) + chr$(13)
   t$ = t$ + "Covariance du couple (x,y) : cov(x,y) = " + str$(int(cov*1E3)/1E3) +chr$(13)
   t$ = t$ + "Coefficient de corrélation : c. corr = " + str$(int(corr*1E3)/1E3) + chr$(13)
   t$ = t$ + "écart-type(x) = " + str$(int(ect_x*1E3)/1E3) +chr$(13)
   t$ = t$ + "écart-type(y) = " + str$(int(ect_y*1E3)/1E3)
   caption 20, t$
END_SUB
rem ============================================================================
FNC Permuter(k)
   dim_local i,q,li,aux
   q = 0 : li = k
   while a(li,k) = 0
       q = q+1 : li = q+k
       if li = n+1 then result 0 : exit_fnc
   end_while
   for i = 1 to n+1
       aux = a(k,i) : a(k,i) = a(li,i) : a(li,i) = aux
   next i
   result a(k,k)
END_FNC
rem ============================================================================
SUB Singulier()
   message "Pas de chance, ce système est singulier ..."
END_SUB
rem ============================================================================
' Tracer la droite y = ax + b
SUB Tracer()
   dim_local a,b
   a = int(a1*1E3)/1E3
   b = int(b1*1E3)/1E3
   2d_pen_color 0,0,255
   2d_line xo-10,yo-npt-b*unite,xo+npt*unite,yo-npt-unite*(a*npt+b)  
END_SUB
rem ============================================================================

chart 1
chart_point 1
item_add 1,"1,0.5"
item_add 1,"2,1.8"
item_add 1,"3,1.8"
item_add 1,"4,2.2"
item_add 1,"5,2.7"
item_add 1,"6,3.3"
item_add 1,"7,3.5"
item_add 1,"8.2,4"
item_add 1,"9,4.7"
item_add 1,"10,5"
item_add 1,"11,6"
item_add 1,"12,6"
item_add 1,"13,7"
item_add 1,"14,7.8"
item_add 1,"15,9"
item_add 1,"16,9.5"

chart_line 1
item_add 1,"0,-0.076"
item_add 1,"17,9.56"

type TModel =
  (REG_LIN,    { y = a + b.x }
   REG_POL2,   { y = a + b.x + c.x^2 }
   REG_POL3,   { y = a + b.x + c.x^2 + d.x^3 }
   REG_FRAC1,  { y = a.x / (1 + b.x) }
   REG_FRAC2,  { y = (a + b.x) / (1 + c.x) }
   REG_EXPO1,  { y = a.exp(-b.x) }
   REG_EXPO2,  { y = a + b.exp(-c.x) }
   REG_LOGIS1, { y = a / [1 + exp(-b.x + c)] }
   REG_LOGIS2, { y = a + (b - a) / [1 + exp(-c.x + d)] }
   REG_POWER1, { y = a.x^b }
   REG_POWER2, { y = a / [1 + (b / x)^c] }
   REG_USER);  { Fonction definie par l'utilisateur }