96
ch−¬ng 15
phÐp néi suy vµ mÞn ho¸ ®-êng cong
Trong c¸c lÜnh vùc øng dông sè, nhiÖm vô cña chóng ta lµ ph¶i biÓu diÔn sè liÖu, th−êng lµ c¸c
sè ®o b»ng c¸c chøc n¨ng ph©n tÝch. Cã hai c¸ch gi¶i quyÕt vÊn ®Ò nµy, trong ph−¬ng ph¸p nèi ®iÓm
(interpolation) th× d÷ liÖu ®−îc coi lµ ®óng vµ c¸i chóng ta cÇn lµ c¸ch biÓu diÔn d÷ liÖu kh«ng n»m
gi÷a c¸c gi¸ trÞ ®o ®−îc, theo ph−¬ng ph¸p thø hai gäi lµ ph−¬ng ph¸p mÞn ho¸ ®õng cong (curve
fitting or regression), b¹n t×m mét ®õng cong kh«ng g·y khóc mµ phï hîp nhÊt víi d÷ liÖu ®· cã, nh-
−ng kh«ng cÇn thiÕt ph¶i ®i qua mét c¸ch chÝnh x¸c bÊt kú mét ®iÓn nµo trªn b¶ng sè liÖu. H15.1
minh ho¹ hai ph−¬ng ph¸p trªn, ch÷ o ®¸nh dÊu c¸c ®iÓm biÓu diÔn d÷ liÖu, c¸c ®o¹n th¼ng b»ng nÐt
liÒn nèi c¸c ®−êng biÓu diÔn d÷ liÖu l¹i víi nhau theo ph¬ng ph¸p nèi ®iÓm cßn ®−êng chÊm chÊm lµ
mét ®õng cong vÏ theo ph−¬ng ph¸p mÞn ho¸ d÷ liÖu.
15.1 MÞn ho¸ ®−êng cong
Ph−¬ng ph¸p mÞn ho¸ ®−êng cong liªn quan ®Õn viÖc tr¶ lêi hai c©u hái c¬ b¶n, ®ã lµ ®−êng cong
thÕ nµo th× phï hîp víi d÷ liÖu nhÊt vµ c©u hái thø hai lµ ph¶i sö dông lo¹i ®−êng cong nµo. “Phï hîp
nhÊt” cã thÓ hiÓu theo nhiÒu c¸ch vµ do ®ã cã nhiÒu ®−êng cong, v× vËy chóng ta ph¶i b¾t ®Çu tõ ®©u?.
NÕu “phï hîp nhÊt” lµ gi¶m nhá ®Õn møc tèi thiÓu tæng sai sè qu©n ph−¬ng t¹i mçi ®iÓm biÓu diÔn d÷
liÖu, so víi gi¸ trÞ t−¬ng øng trªn ®−êng cong th× ®−êng cong phï hîp nhÊt sÏ lµ mét ®−êng th¼ng vÒ
mÆt to¸n mµ nãi ph−¬ng ph¸p nµy ®−îc gäi lµ ph−¬ng ph¸p xÊp xØ ®a thøc. NÕu nh− kh¸i niÖm nµy
cßn khã hiÓu ®èi víi b¹n th× xin h·y xem l¹i h×nh 15.1 kho¶ng c¸ch theo chiÒu däc gi÷a ®−êng cong
d÷ liÖu vµ c¸c ®iÓm biÓu diÔn d÷ liÖu gäi lµ sai sè cña ®iÓm ®ã, b×nh ph−¬ng kho¶ng c¸ch nµy lªn vµ
céng tÊt c¶ chóng l¹i ta ®îc tæng b×nh ph−¬ng sai sè. §−êng cong chÊm chÊm lµ ®−êng cong lµm cho
b×nh ph−¬ng sai sè lµ nhá nhÊt vµ ®−îc gäi lµ ®−êng cong phï hîp nhÊt. Tõ “qu©n ph−¬ng bÐ nhÊt” lµ
c¸ch nãi t¾t cña côm tõ “Tæng b×nh ph¬ng sai sè bÐ nhÊt”.
H×nh 15.1
97
Trong MATLAB hµm
polyfit
sÏ gi¶i quyÕt vÊn ®Ò xÊp xØ ®−êng cong qu©n ph−¬ng bÐ nhÊt. §Ó
minh ho¹ cho viÖc sö dông hµm nµy, chóng ta h·y b¾t ®Çu b»ng c¸c d÷ liÖu ®· cã ë trong h×nh vÏ.
>> x = [0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1];
>> y =[-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2];
§Ó sö dông hµm
polyfit
, chóng ta ph¶i truyÒn cho nã d÷ liÖu trªn vµ bËc cña ®a thøc mµ chóng ta
muèn phï hîp víi d÷ liÖu, nÕu chóng ta chän bËc n lµ 1 th× ®−êng cong xÊp xØ gÇn nhÊt sÏ lµ ®−êng
th¼ng. Ph−¬ng ph¸p nµy ®−îc gäi lµ ph−¬ng ph¸p xÊp xØ tuyÕn tÝnh. MÆt kh¸c nÕu chóng ta chon n=2
th× chóng ta sÏ t×m ®−îc mét tam thøc bËc hai. VÝ dô:
>> n = 2;
>> p = polyfit(x,y,n)
p =
-9.8108 20.1293 -0.0317
KÕt qu¶ cña
polyfit
lµ mét vector biÓu diÔn hÖ sè cña mét ®a thøc bËc hai. ë ®©y ®a thøc ®ã lµ
y= -9.8108x2+20.1293x-0.0317. §Ó so s¸nh møc ®é xÊp xØ cña ®a thøc víi c¸c ®iÓm d÷
liÖu chóng ta h·y vÏ hai ®−êng:
>> xi = linspace(0,1,100);
Dßng nµy ®Ó t¹o ra d÷ liÖu trôc x ®Ó chuÈn bÞ vÏ ®a thøc
>> z = polyval(p,xi)
Dßng nµy gäi hµm
polyval
cña MATLAB ®Ó tÝnh gi¸ trÞ cña ®a thøc p t¹i c¸c ®iÓm xi
>> plot(x,y,'-o',xi,z,':')
VÏ c¸c ®iÓm cã to¹ ®é lµ x vµ y, ®¸nh dÊu c¸c ®iÓm nµy b»ng ch÷ ‘o’ sau ®ã nèi c¸c ®iÓm nµy b»ng
c¸c ®o¹n th¼ng. Ngoµi ra nã cßn vÏ d÷ liÖu cña ®a thøc xi vµ z dïng ®−êng chÊm chÊm.
>> xlabel('x'),ylabel('y=f(x)')
>> title('Second Oder Curver Fitting')
T¹o nh·n cho ®−êng cong võa vÏ. KÕt qu¶ cña c¸c lÖnh trªn ®©y lµ mét ®å thÞ ®· ®−îc giíi thiÖu ë
trªn.
ViÖc chän bËc cña ®a thøc kh«ng ph¶i lµ ngÉu nhiªn, nÕu cã hai ®iÓm th× x¸c ®Þnh mét ®−êng
th¼ng, tøc lµ mét ®a thøc bËc nhÊt, ba ®iÓm th× x¸c ®Þnh mét parabol bËc hai. Cø nh− vËy, ®Ó x¸c ®Þnh
mét ®−êng cong bËc n, cÇn cã n+1 ®iÓm. V× vËy, ë trong vÝ dô tr−íc cã 11 ®iÓm d÷ liÖu, chóng ta cã
thÓ chän bËc cña ®a thøc lµ tõ 1 ®Õn 10. Tuy nhiªn, do tÝnh chÊt sè häc cña c¸c ®a thøc bËc cao rÊt
phøc t¹p nªn b¹n kh«ng nªn chän bËc cña ®a thøc lín h¬n møc cÇn thiÕt. Ngoµi ra khi bËc cña ®a thøc
t¨ng lªn th× sù xÊp xØ cµng kÐm h¬n, v× vËy c¸c ®a thøc bËc cao cã thÓ bÞ ®¹o hµm nhiÒu lÇn tr−íc khi
®¹o hµm cña chóng b»ng kh«ng. VÝ dô cho mét ®a thøc bËc 10:
>> pp = polyfit(x,y,10)
pp =
1.0e+006 *
Columns 1 through 7
-0.4644 2.2965 -4.8773 5.8233 -4.2948 2.0211 -
0.6032
Columns 8 through 11
98
0.1090 -0.0106 0.0004 -0.0000
>> format short e % change display format
>> pp.' % display polynomial coefficients as a column
ans =
-4.6436e+005
2.2965e+006
-4.8773e+006
5.8233e+006
-4.2948e+006
2.0211e+006
-6.0322e+005
1.0896e+005
-1.0626e+004
4.3599e+002
-4.4700e-001
L−u ý kÝch th−íc cña vector hÖ sè ®a thøc trong tr−êng hîp nµy so víi ®−êng cong bËc hai tr−íc
®©y, ®ång thêi còng lu ý sù kh¸c nhau gi÷a sè h¹ng nhá nhÊt vµ sè h¹ng lín nhÊt trong ®a thøc vµo
kho¶ng 107. H·y thö vÏ ®−êng cong nµy vµ so s¸nh víi d÷ liÖu gèc vµ víi ®−êng cong bËc hai.
>> zz = polyval(pp,xi); % evalute 10th order polynomial
>> plot(x,y,'o',xi,z,’:’,xi,zz) % plot data
>> xlabel('x'),ylabel('y=f(x)')
>> title('2nd and 10th Order Curver Fitting')
H×nh 15.2
99
Trªn h×nh 15.2, d÷ liÖu gèc ®−îc ®¸nh dÊu o, ®−êng cong bËc hai ®−îc vÏ b»ng nÐt chÊm chÊm,
cßn ®−êng cong bËc 10 ®−îc vÏ b»ng nÐt ®Ëm. §Ó ý ®Õn nÐt gîn sãng xuÊt hiÖn gi÷a c¸c ®iÓm d÷ liÖu
bªn phÝa tr¸i vµ bªn phÝa ph¶i cña ®−êng cong bËc 10. Dùa vµo ®å thÞ nµy th× râ rµng r»ng c¸i chiÕt lý
cµng nhiÒu cµng tèt kh«ng thÓ ¸p dông ®−îc ë ®©y.
15.2 Nèi ®iÓm mét chiÒu
Nh− ®· giíi thiÖu th× nèi ®iÓm ®−îc ®Þnh nghÜa nh− lµ mét ph−¬ng ph¸p dù ®o¸n gi¸ trÞ cña hµm
gi÷a nh÷ng ®iÓm cho tr−íc. Nèi ®iÓm lµ mét c«ng cô h÷u hiÖu khi chóng ta kh«ng thÓ nhanh chãng
tiÝnh ®−îc gi¸ trÞ cña hµm t¹i c¸c ®iÓm trung gian. Ph−¬ng ph¸p nµy ®−îc sö dông réng r·i ®èi víi d÷
liÖu lµ gi¸ trÞ cña c¸c phÐp ®o thùc nghiÖm hoÆc lµ kÕt qu¶ cña c¸c chuçi tÝnh to¸n dµi. Cã thÓ vÝ dô
®¬n gi¶n nhÊt cña viÖc nèi ®iÓm chÝnh lµ ph−¬ng ph¸p vÏ tõng ®iÓm cña MATLAB, tøc lµ vÏ nh÷ng
®o¹n th¼ng nèi nh÷ng ®iÓm d÷ liÖu liªn tiÕp ®Ó t¹o lªn mét ®å thÞ.
§©y lµ ph−¬ng ph¸p nèi ®iÓm tuyÕn tÝnh, nã cho r»ng c¸c gi¸ trÞ cña hµm n»m gi÷a hai ®iÓm cho
tr−íc sÏ r¬i vµo kho¶ng gi÷a hai ®Çu cña ®o¹n th¼ng nèi hai ®iÓm ®ã. HiÓn nhiªn lµ khi sè l−îng c¸c
®iÓm d÷ liÖu t¨ng lªn vµ kho¶ng c¸ch gi÷a chóng gi¶m ®i th× ph−¬ng ph¸p nèi ®iÓm tuyÕn tÝnh cµng
trë lªn chÝnh x¸c.
>> x1 = linspace(0,2*pi,60);
>> x2 = linspace(0,2*pi,6);
>> plot(x1,sin(x1),x2,sin(x2),'-')
>> xlabel('x'),ylabel('sin(x)')
>> title('Linear Interpolation')
H×nh 15.3
C¶ hai ®å thÞ cïng vÏ mét hµm sine nh−ng ®å thÞ 60 ®iÓm th× mÞn h¬n ®å thÞ 6 ®iÓm.
100
Còng gièng nh− ph−¬ng ph¸p xÊp xØ ho¸ ®−êng cong, ë ®©y chóng ta còng ph¶i thùc hiÖn mét sè
lùa chän, cã rÊt nhiÒu c¸ch ®Ó nèi hai ®iÓm, tuú thuéc vµo gi¶ ®Þnh mµ chóng ta ®· lùa chän. H¬n n÷a
chóng ta cã thÓ nèi c¸c ®iÓm trong kh«ng gian kh«ng ph¶i lµ mét chiÒu. Nãi nh− thÕ nÕu b¹n cã d÷
liÖu ph¶n ¸nh mét hµm phô thuéc vµo hai biÕn z=f(x,y), b¹n cã thÓ nèi gi¸ trÞ n»m gi÷a hai ®iÓm cã x
vµ y kh¸c nhau ®Ó t×m ra gi¸ trÞ trung gian cña hai ®iÓm. MATLAB cung cÊp mét sè hµm ®Ó nèi lµ :
interp1
nèi c¸c d÷ liÖu mét chiÒu,
interp2
nèi c¸c d÷ liÖu hai chiÒu,
interp3
nèi c¸c d÷ liÖu ba chiÒu,
interpn
nèi c¸c d÷ liÖu cã sè chiÒu lín h¬n 3.
Sau ®©y chóng ta sÏ xem xÐt c¸c d÷ liÖu mét vµ hai chiÒu. §Ó minh ho¹ viÖc nèi d÷ liÖu mét chiÒu,
h·y xÐt vÝ dô sau, kh¶ n¨ng cña thÝnh gi¸c, vÝ dô nh møc ©m thanh bÐ nhÊt hay cßn gäi lµ ng−ìng
nghe cña tai ngêi thay ®æi theo tÇn sè, d÷ liÖu do ng−êi thèng kª ®−îc cho nh− sau:
>> Hz = [20:10:100 200:100:1000 1500 2000:1000:10000];
>> % Frequencies in Hertz
>> spl = [76 66 59 54 49 46 43 40 38 22 ...
14 9 6 3.5 2.5 1.4 0.7 0 -1 -3 ...
-8 -7 -2 2 7 9 11 12];
>> % sound pressure level in dB
Ng−ìng nghe ®−îc chuÈn ho¸ b»ng 0dB t¹i tÇn sè 1000Hz, bëi v× tÇn sè tr¶i trong mét d¶i rÊt réng
nªn khi vÏ c¸c ®iÓm d÷ liÖu chóng ta logarithm ho¸ trôc x.
>> semilogx(Hz,spl,'-o')
>> xlabel('Frequency, Hz')
>> ylabel('Relative Sound Presure Level1, dB')
>> title('Threshold of Human Hearing')
Dùa vµo h×nh 15.4 ta thÊy tai ng−êi nh¹y c¶m hÕt ®èi víi c¸c ©m thanh trong kho¶ng 3kHz. Dùa
vµo c¸c sè liÖu nµy, chóng ta h·y dù ®o¸n ng−ìng nghe ë tÇn sè 2,5kHz b»ng mét vµi c¸ch kh¸c nhau.
>> s = interp1(Hz,spl,2.5e3) %linear interpolation
s =
-5.5000e+000
>> s = interp1(Hz,spl,2.5e3,'linear') %linear interpolation again
s =
-5.5000e+000
>> s = interp1(Hz,spl,2.5e3,'cubic') % cubic interpolation
s =
-5.8690e+000
>> s = interp1(Hz,spl,2.5e3,'spline') % spline interpolation
s =
-5.8690e+000
>> s = interp1(Hz,spl,2.5e3,'nearest')% nearest-neighbor
s =
-8
H·y ®Ó ý ®Õn sù kh¸c nhau trong c¸c kÕt qu¶, hai gi¸ trÞ ®Çu tiªn tr¶ vÒ mét c¸ch chÝnh x¸c gi¸ trÞ
®−îc vÏ ë trªn h×nh t¹i tÇn sè 2,5kHz bëi v× MATLAB ®· nèi c¸c ®iÓm mét c¸ch tuyÕn tÝnh gi÷a c¸c
®iÓm d÷ liÖu trªn ®å thÞ c¸c ®−êng cong ®a thøc, vÝ dô nh− ®a thøc bËc 3 sÏ xÊp xØ ho¸ c¸c ®iÓm trªn
®å thÞ theo c¸c c¸ch kh¸c nhau, kÕt qu¶ lµ c¸c ®−êng cong nµy t−¬ng ®èi phï hîp víi c¸c d÷ liÖu mµ
nã ®i qua trªn ®å thÞ nh−ng kh¸c biÖt kh¸ xa so víi ph−¬ng ph¸p nèi b»ng ®−êng th¼ng.