Estimation moindres carrés - validation

Version imprimable

Bonjour,

J'ai effectué une estimation de moindres carrés pour obtenir les constantes a0, a1, a2, b1, b2,c1 de la fonction suivante:
FF=a0 + a1 * altitude+ a2 * altitude**2+ b1 * FN + b2 * FN+ c1 * altitude * FN
J'avais pour base des valeurs avec un pas de 100.

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 
altitude = [28000.0, 28100.0, 28200.0, 28300.0, 28400.0, 28500.0, 28600.0, 28700.0, 28800.0, 28900.0, 29000.0, 29100.0, 29200.0, 29300.0, 29400.0, 29500.0, 29600.0, 29700.0, 29800.0, 29900.0, 30000.0, 30100.0, 30200.0, 30300.0, 30400.0, 30500.0, 30600.0, 30700.0, 30800.0, 30900.0, 31000.0, 31100.0, 31200.0, 31300.0, 31400.0, 31500.0, 31600.0, 31700.0, 31800.0, 31900.0, 32000.0, 32100.0, 32200.0, 32300.0, 32400.0, 32500.0, 32600.0, 32700.0, 32800.0, 32900.0, 33000.0, 33100.0, 33200.0, 33300.0, 33400.0, 33500.0, 33600.0, 33700.0, 33800.0, 33900.0, 34000.0, 34100.0, 34200.0, 34300.0, 34400.0, 34500.0, 34600.0, 34700.0, 34800.0, 34900.0, 35000.0, 35100.0, 35200.0, 35300.0, 35400.0, 35500.0, 35600.0, 35700.0, 35800.0, 35900.0, 36000.0, 36100.0, 36200.0, 36300.0, 36400.0, 36500.0, 36600.0, 36700.0, 36800.0, 36900.0, 37000.0, 37100.0, 37200.0, 37300.0, 37400.0, 37500.0, 37600.0, 37700.0, 37800.0, 37900.0, 38000.0, 38100.0, 38200.0, 38300.0, 38400.0, 38500.0, 38600.0, 38700.0, 38800.0, 38900.0, 39000.0]
 
FN = [2965.0, 2955.0, 2946.0, 2937.0, 2927.0, 2918.0, 2908.0, 2899.0, 2890.0, 2880.0, 2871.0, 2861.0, 2852.0, 2842.0, 2831.0, 2820.0, 2810.0, 2799.0, 2788.0, 2778.0, 2767.0, 2757.0, 2747.0, 2737.0, 2727.0, 2717.0, 2707.0, 2697.0, 2687.0, 2677.0, 2667.0, 2657.0, 2647.0, 2637.0, 2627.0, 2617.0, 2608.0, 2598.0, 2588.0, 2578.0, 2568.0, 2559.0, 2549.0, 2539.0, 2529.0, 2520.0, 2510.0, 2500.0, 2491.0, 2481.0, 2471.0, 2462.0, 2452.0, 2442.0, 2433.0, 2423.0, 2414.0, 2404.0, 2395.0, 2385.0, 2376.0, 2366.0, 2356.0, 2346.0, 2336.0, 2325.0, 2315.0, 2305.0, 2295.0, 2285.0, 2275.0, 2265.0, 2255.0, 2245.0, 2235.0, 2225.0, 2215.0, 2205.0, 2196.0, 2186.0, 2176.0, 2165.0, 2155.0, 2144.0, 2134.0, 2123.0, 2113.0, 2103.0, 2092.0, 2082.0, 2072.0, 2062.0, 2052.0, 2042.0, 2032.0, 2022.0, 2012.0, 2002.0, 1992.0, 1983.0, 1973.0, 1963.0, 1954.0, 1944.0, 1935.0, 1925.0, 1916.0, 1906.0, 1897.0, 1888.0, 1878.0]
 
fuelflow = [2022.13, 2012.3550000000002, 2006.226, 2000.0970000000002, 1993.287, 1987.1580000000001, 1980.3480000000002, 1971.3200000000002, 1965.2, 1958.4, 1952.2800000000002, 1945.4800000000002, 1939.3600000000001, 1932.5600000000002, 1925.0800000000002, 1917.6000000000001, 1907.9900000000002, 1900.5210000000002, 1893.0520000000001, 1886.2620000000002, 1878.7930000000001, 1872.0030000000002, 1865.2130000000002, 1858.4230000000002, 1851.633, 1844.843, 1838.053, 1831.2630000000001, 1824.4730000000002, 1817.6830000000002, 1810.893, 1804.103, 1797.313, 1790.5230000000001, 1783.7330000000002, 1776.9430000000002, 1770.832, 1764.0420000000001, 1757.2520000000002, 1750.4620000000002, 1743.672, 1735.0020000000002, 1728.2220000000002, 1721.442, 1714.662, 1708.5600000000002, 1701.7800000000002, 1695.0000000000002, 1688.8980000000001, 1682.1180000000002, 1675.3380000000002, 1669.236, 1662.4560000000001, 1655.6760000000002, 1649.574, 1642.794, 1634.278, 1627.508, 1621.4150000000002, 1614.6450000000002, 1608.5520000000001, 1601.7820000000002, 1595.0120000000002, 1588.2420000000002, 1581.4720000000002, 1574.025, 1567.255, 1560.4850000000001, 1553.7150000000001, 1546.9450000000002, 1540.1750000000002, 1533.4050000000002, 1526.635, 1519.865, 1513.095, 1504.1000000000001, 1497.3400000000001, 1490.5800000000002, 1484.496, 1477.736, 1470.976, 1463.5400000000002, 1456.7800000000002, 1449.344, 1442.584, 1435.1480000000001, 1428.3880000000001, 1421.6280000000002, 1414.192, 1407.432, 1400.672, 1393.912, 1387.152, 1382.4340000000002, 1375.664, 1368.894, 1362.124, 1355.354, 1348.584, 1342.491, 1335.721, 1328.951, 1322.8580000000002, 1316.0880000000002, 1309.9950000000001, 1303.2250000000001, 1297.132, 1290.362, 1284.269, 1278.1760000000002, 1271.4060000000002]
 
def func(alt, fn, p):
    a0, a1, a2, b1, b2, c1 = p
    return a0 + a1 * alt + a2 * (alt)** 2 + b1 * fn + b2 * (fn)**2 + c1 * alt  * fn
 
def residuals(p):
    return [ff - func(a, fn, p) for ff, a, fn in zip(fuelflow, altitude, FN)]
 
plsq, ler = leastsq(residuals, x0=[1] * 6)
print('plsq : \n', plsq)
 
FuelFlow_estime = [func(a, fn, plsq) for a, fn in zip(altitude, FN)]
print('FuelFlow_estime : \n', FuelFlow_estime)

J'obtiens les coefficients suivants:

Code:

1
2
3
 
[  1.10906252e+06,  -3.83274171e+01,   3.31154773e-04,
        -3.84702938e+02,   3.34848674e-02,   6.65828681e-03]

Quand ensuite j'utilise ces coefficients pour calculer des FF j'ai des valeurs avec un ordre de grandeur tout à fait incohérent.
Par exemple, avec le code suivant:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
 
def Calcul_FF (m,FL):
    #on divise le Fn par 10 car la regression des moindres carres avait été effectuée avec des valeurs en daN#
    Fn = 5067.898977626907
    a0 = 1.10906252e+06
    a1 = -3.83274171e+01
    a2 = 3.31154773e-04
    b1 = -3.84702938e+02
    b2 = 3.34848674e-02
    c1 = 6.65828681e-03
    # le fuel flow est donné en kg/h #
    FF = a0 + a1 * FL * 100 + a2 *  (FL*100)**2 + b1 * Fn + b2 * (Fn**2) + c1 * (FL * 100 * Fn)
    return FF
 
print ("FF =",Calcul_FF(83000,360))

j'obtiens :

Code:

1
2
 
FF = 283595.30298786773

là où je devrais retrouver un ordre de grandeur de 10^3 max..

Je voudrais donc ploter l'estimation et la courbe que je calcule moi pour voir si elle fait "des dents de scie" autour de l'estimation. Et j'aimerais calculer le coefficient de corrélation de l'estimation: comment faire ?

merci pour votre aide !

Salut,

Intéressant. Quel est l'intérêt d'avoir une function func() et une autre function Calcul_FF() ? Ne sont-elles pas toutes deux en charge de calculer la meme chose? Ton fitting m'a l'air bon vu comme ça. En faisant un test (en copiant ton code), FuelFlow_estime donne de bonnes valeurs.

J'ai pas regardé dans le detail tout ce que tu as fait, mais l'erreur peut venir de ta conversion d'unité. Fais la conversion sur FF et non sur Fn dans ta function Cacul_FF().

Sinon, pour ce qui est du coefficient de correlation, tu peux utiliser la function proposée par numpy, numpy.corrcoeff(x,y). Des details sur son utilisation se trouvent ici:

http://docs.scipy.org/doc/numpy-1.6.....corrcoef.html

Voici ton code en rajoutant le calcul du r²:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
 
from __future__ import division
from scipy.optimize import leastsq
import numpy as np
 
# Data
altitude = [28000.0, 28100.0, 28200.0, 28300.0, 28400.0, 28500.0, 28600.0, 28700.0, 28800.0, 28900.0, 29000.0, 29100.0, 29200.0, 29300.0, 29400.0, 29500.0, 29600.0, 29700.0, 29800.0, 29900.0, 30000.0, 30100.0, 30200.0, 30300.0, 30400.0, 30500.0, 30600.0, 30700.0, 30800.0, 30900.0, 31000.0, 31100.0, 31200.0, 31300.0, 31400.0, 31500.0, 31600.0, 31700.0, 31800.0, 31900.0, 32000.0, 32100.0, 32200.0, 32300.0, 32400.0, 32500.0, 32600.0, 32700.0, 32800.0, 32900.0, 33000.0, 33100.0, 33200.0, 33300.0, 33400.0, 33500.0, 33600.0, 33700.0, 33800.0, 33900.0, 34000.0, 34100.0, 34200.0, 34300.0, 34400.0, 34500.0, 34600.0, 34700.0, 34800.0, 34900.0, 35000.0, 35100.0, 35200.0, 35300.0, 35400.0, 35500.0, 35600.0, 35700.0, 35800.0, 35900.0, 36000.0, 36100.0, 36200.0, 36300.0, 36400.0, 36500.0, 36600.0, 36700.0, 36800.0, 36900.0, 37000.0, 37100.0, 37200.0, 37300.0, 37400.0, 37500.0, 37600.0, 37700.0, 37800.0, 37900.0, 38000.0, 38100.0, 38200.0, 38300.0, 38400.0, 38500.0, 38600.0, 38700.0, 38800.0, 38900.0, 39000.0]
FN = [2965.0, 2955.0, 2946.0, 2937.0, 2927.0, 2918.0, 2908.0, 2899.0, 2890.0, 2880.0, 2871.0, 2861.0, 2852.0, 2842.0, 2831.0, 2820.0, 2810.0, 2799.0, 2788.0, 2778.0, 2767.0, 2757.0, 2747.0, 2737.0, 2727.0, 2717.0, 2707.0, 2697.0, 2687.0, 2677.0, 2667.0, 2657.0, 2647.0, 2637.0, 2627.0, 2617.0, 2608.0, 2598.0, 2588.0, 2578.0, 2568.0, 2559.0, 2549.0, 2539.0, 2529.0, 2520.0, 2510.0, 2500.0, 2491.0, 2481.0, 2471.0, 2462.0, 2452.0, 2442.0, 2433.0, 2423.0, 2414.0, 2404.0, 2395.0, 2385.0, 2376.0, 2366.0, 2356.0, 2346.0, 2336.0, 2325.0, 2315.0, 2305.0, 2295.0, 2285.0, 2275.0, 2265.0, 2255.0, 2245.0, 2235.0, 2225.0, 2215.0, 2205.0, 2196.0, 2186.0, 2176.0, 2165.0, 2155.0, 2144.0, 2134.0, 2123.0, 2113.0, 2103.0, 2092.0, 2082.0, 2072.0, 2062.0, 2052.0, 2042.0, 2032.0, 2022.0, 2012.0, 2002.0, 1992.0, 1983.0, 1973.0, 1963.0, 1954.0, 1944.0, 1935.0, 1925.0, 1916.0, 1906.0, 1897.0, 1888.0, 1878.0]
fuelflow = [2022.13, 2012.3550000000002, 2006.226, 2000.0970000000002, 1993.287, 1987.1580000000001, 1980.3480000000002, 1971.3200000000002, 1965.2, 1958.4, 1952.2800000000002, 1945.4800000000002, 1939.3600000000001, 1932.5600000000002, 1925.0800000000002, 1917.6000000000001, 1907.9900000000002, 1900.5210000000002, 1893.0520000000001, 1886.2620000000002, 1878.7930000000001, 1872.0030000000002, 1865.2130000000002, 1858.4230000000002, 1851.633, 1844.843, 1838.053, 1831.2630000000001, 1824.4730000000002, 1817.6830000000002, 1810.893, 1804.103, 1797.313, 1790.5230000000001, 1783.7330000000002, 1776.9430000000002, 1770.832, 1764.0420000000001, 1757.2520000000002, 1750.4620000000002, 1743.672, 1735.0020000000002, 1728.2220000000002, 1721.442, 1714.662, 1708.5600000000002, 1701.7800000000002, 1695.0000000000002, 1688.8980000000001, 1682.1180000000002, 1675.3380000000002, 1669.236, 1662.4560000000001, 1655.6760000000002, 1649.574, 1642.794, 1634.278, 1627.508, 1621.4150000000002, 1614.6450000000002, 1608.5520000000001, 1601.7820000000002, 1595.0120000000002, 1588.2420000000002, 1581.4720000000002, 1574.025, 1567.255, 1560.4850000000001, 1553.7150000000001, 1546.9450000000002, 1540.1750000000002, 1533.4050000000002, 1526.635, 1519.865, 1513.095, 1504.1000000000001, 1497.3400000000001, 1490.5800000000002, 1484.496, 1477.736, 1470.976, 1463.5400000000002, 1456.7800000000002, 1449.344, 1442.584, 1435.1480000000001, 1428.3880000000001, 1421.6280000000002, 1414.192, 1407.432, 1400.672, 1393.912, 1387.152, 1382.4340000000002, 1375.664, 1368.894, 1362.124, 1355.354, 1348.584, 1342.491, 1335.721, 1328.951, 1322.8580000000002, 1316.0880000000002, 1309.9950000000001, 1303.2250000000001, 1297.132, 1290.362, 1284.269, 1278.1760000000002, 1271.4060000000002]
 
def func(alt, fn, p):
    a0, a1, a2, b1, b2, c1 = p
    return a0 + a1 * alt + a2 * (alt)** 2 + b1 * fn + b2 * (fn)**2 + c1 * alt  * fn
 
def residuals(p):
    return [ff - func(a, fn, p) for ff, a, fn in zip(fuelflow, altitude, FN)]
 
plsq, ler = leastsq(residuals, x0=[1] * 6)
print('plsq : \n', plsq)
 
FuelFlow_estime = [func(a, fn, plsq) for a, fn in zip(altitude, FN)]
print('FuelFlow_estime : \n', FuelFlow_estime)
 
print('Correlation Coeff: \n', np.corrcoef(fuelflow, FuelFlow_estime)[0,1]**2)

Sinon pour ce qui est de visualizer, je te conseil de tracer les points de fuelflow et une courbe faite avec FuelFlow_estime.

Tu peux aller voir cette page qui peut te donner des idées:
http://nbviewer.ipython.org/github/d...eFitting.ipynb

Ciao

1 pièce(s) jointe(s)

Salut.
Je pense aussi qu'il y a un problème de conversion d'unités. Ça n'est pas très clair pour qu'on puisse t'aider plus.

La fonction scipy.optimize.curve_fit est mieux adaptée, et plus facile à utiliser que leastsq. Elle retourne une estimation de la covariante. (On peut l'avoir avec leastsq, avec l'argument full_output=True, mais il faire un peu plus de travail).

J'utilise matplotlib pour visualiser. Par exemple :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt ... fig = plt.figure() ax = fig.add_subplot(111, projection='3d') for c, m, d in (('b', 'o', FuelFlow_estime), ('r', '^', fuelflow)): for coords in zip(altitude, FN, d): ax.scatter(*coords, c=c, marker=m) x = np.linspace(27000, 40000, 20) y = np.linspace(1800, 3000, 20) X, Y = np.meshgrid(x, y) Z = func(X, Y, plsq) ax.plot_wireframe(X, Y, Z,) ax.set_xlabel('altitude') ax.set_ylabel('FN') ax.set_zlabel('FuelFlow') plt.show()
C'est pas super pour faire un fit d'avoir les points sur la même droite.
Pièce jointe 147994

09/06/2014, 13h10
Julien N

Salut,

Citation:

Envoyé par __dardanos__

Je pense aussi qu'il y a un problème de conversion d'unités. Ça n'est pas très clair pour qu'on puisse t'aider plus.

C'est même certain. Dans la formule FF=a0 + a1 * altitude+ a2 * altitude**2+ b1 * FN + b2 * FN+ c1 * altitude * FN, les coefficients a0, a1, a2, b1, b2 et c1 ne sont pas adimensionnels.

Joli graph !

Ju