?? sm6.htm
字號:
<html><head><title>幫助</title></head>
<body>
<center>
<table width=780 border=0 cellpadding=3 cellspacing=0>
<tr><td style='line-height:150%;font-size:16px;text-indent: 20'>
<p><a href="readme.htm">返回</a></p>
<p align=center><b>第六章：離散序列的直線擬合算法</b></p>
<p align=center>許劍偉 2008-12-2 于莆田十中</p>
<p><b>問題提出：</b>已知古代朔日表，如何找到一條直線 y(x) = k·x + b 擬合朔日表，然后用[y(x)]表示朔日，x表示朔日序數(x=0,1,2……)，k和b是待求系數。</p>
<p><b>數學模型：</b>設某一組實測的離散數據序列D(x)，其中x = 0,1,2……。現求解一條直線 y(x) = k·x + b，使y(x)在 x = 0,1,2……的值與實測值的誤差大等于0小于a。對于上述問題，a取值為1，遺憾的是這個問題不能使用最小二乘法直接求解，我們考慮其它方法。</p>
<p><b>數學模型的圖解描述</b>：</p>
<p align=center><img src=img/a001.gif></p>
<p>圖中D(x)是已知的離散序列，E(x)= D(x)+ a，擬合直線為y = k·x + b，找出適當的k和b使得直線y(x)完全穿過紅色區域。</p>
<p><b>解題思路：</b></p>
<p>（1）y(x)是未知的，先取個k和b的估值，如圖中示意，可取k=1，b=0做為初始估值。</p>
<p>（2）找b的值：平移直線（改變b值），直到直線全部大等于D(x)，設b值為b1時剛好超過所有的點，其中臨界超過的點是D(x1)。繼續平移，找到直線剛好全部小于E(x)時b的取值，記為b2。</p>
<p>如果b2>b1，顯然有解，若k值保持不變，b取值（b1+b2）/ 2，那么b允許的誤差為h = （b2-b1）/2，查找過程結束。如果b2<=b1則繼續以下步聚。</p>
<p>（3）找k的值：以點D(x1)為固定點旋轉直線，找出直線不超過E(x)時k的取值范圍（k1，k2）。如果k2<=k1則無解，查找過程結束。否則，k取（k1+k2）/ 2。</p>
<p>（4）從第2步開始重復以上過程，直到成功為止。重復的過程中，k的取值是在原有基礎上進一步縮小范圍的。比如，第一次得到k的范圍是1到2之間，第二次為1.5到3之間，那么k的取值不是取前者也不是取后者，而是取二者的交集（1.5，2），否則以上過程中k的取值無法收斂到有效范圍之內，造成死循環而找不到解。其實，以D(x)中的任意一點為固定點，都可以得到k的上下界值k1和k2，對于一條能夠通過紅區的直線，必然能夠滿足所有k1和k2的約束。</p>
<p>以上得取了b允許的誤差值h，如果b不變，k允許的誤差為b/N，式中N為D(x)序列的總個數。</p>
<p><b>算法實現：</b></p>
<p>設f(x) = D(x) - y(x)，設K(n，m)表示點E(n)到點D(m)連線的斜率。</p>
<p>那么f(x)表示y(x)距D(x)的縱向距離，f(x)+a表示y(x)距E(x)縱向距離。</p>
<p>（1）先取個k和b的初始估值。置k可能的取值范圍為（k1,k2），這個范圍要足夠大。</p>
<p>（2）在x=0，1，2……上遍歷f(x)，求得f(x)的最大值f1和最小值f2及相應的x值x1、x2。那么f2+a就是f(x)+ a的最小值，也就是說f2+a是y(x)到E(x)的最小縱向距離。</p>
<p>（3）如果f1<f2+a則有解，b = b + (f2+a+f1)/2，h = (f2+a-f1)/2，并輸出結果，程序結束。</p>
<p>（4）在x>x1上遍歷K(x，x1)，取得最小值k1（該最小必須須比原來的k1還要小，否則不算）</p>
<p> 在x<x1上遍歷K(x，x1)，取得最大值k2（該最大必值須比原來的k2還要大，否則不算）</p>
<p> 如果k1>=k2則無解，輸出無解報告，程序結束。否則取k = (k1+k2)/2。</p>
<p>（5）從第2步開始重復計算，直到有解。注意， k1和k2是全過程中的最值，而不是某一次遍歷過程中的最值。通常1到6次就能得到結果，如果重復計算10次仍無解，程序也應結束。</p>
<p><b>優化算法：</b></p>
<p>有時我們不走運，程序可以得到滿足條件的直線，但h的值非常小。比如得到h = 0.000002，設N = 10000，那么k允許的誤差為0.0000000002，也就是說k必須取10位有效小數才可保證結果的正確性。所以有必要提高h的值。以下描述本筆者提高h值的方法。</p>
<p>計算前，適當減小a，使得上圖的紅色區域變小，然后出k、b、h，如果k、b有解，那么相應的y(x)定能通過原來較大的紅色區域，這樣b的取值范圍就變得寬裕了，h自然就比較大了。設a的原值記為A，那么只需將上述算法第3步改為：“如果f1<f2+a則有解，b = b + (f2+A+f1)/2，h = (f2+A-f1)/ 2，并輸出結果，程序結束”，等價于“如果f1<f2+a則有解，h = (f2+A-f1)/2，b = b+f1+h，并輸出結果，程序結束”</p>
<p>當然，a的值也不能減小過多。因為某些序列對擬合直線要求非常苛刻，理論上可找到的h值本身就非常小，這是a減小過多可能造成無解。壽星萬年歷的輔助程序中，減小量取值0.0002。</p>
<p><b>范例程序：</b></p>
<hr style='height=1;color:#FF0000'>
<p align=center><b><i>擬合之前先設定參數，以便生成序列。</i></b><br>
D(x)序列生成參數：
k=<input type=text id=Ck value='29.5306' size=10>
b=<input type=text id=Cb value='0.35789' size=10>
N=<input type=text id=Cn value='3000' size=10><br>
擬合前的初使估值：
k=<input type=text id=Ck0 value='29.5' size=10>
b=<input type=text id=Cb0 value='0' size=10>
a=<input type=text id=Ca0 value='0.99998' size=10><br><br>
結果：<input type=text id=Cjg value='' size=60>
<input type=button onclick='test()' value='開始擬合'>
</p>
<hr style='height=1;color:#FF0000'>
<p id=progText></p>
<script language=javascript>
function nihe(r,k,b,a){
//入口參數：擬合計算,k、b為初始估值,r為序列數組,a為上下線的縱向差
//調試時可把a值放寬以便發現問題，取值a=1+1e-7，可靠輸出則取a=1-2e-4。
var i,kk,v, A=1;
var k1=k-1, k2=k+1, f1, f2, x1;
for(kk=0;kk<10;kk++){ //查找次數
f1=-1e9, f2=1e9, x1=-1;
for(i=0;i<r.length;i++){
v=r[i]-(k*i+b); //直線至少要平移的量
if(v>f1) f1=v,x1=i;
if(v<f2) f2=v;
}
if(f2+a-f1>0){ r.gh = (f2+A-f1)/2; b+=f1+r.gh; r.gk = k; r.gb = b; return "ok"; }
for(i=0;i<r.length;i++){ //重新調整k的值
if(i==x1) continue;
v=( r[i]+a - r[x1] ) / (i-x1);
if(i<x1 && v>k1) k1=v;
if(i>x1 && v<k2) k2=v;
if(k1>k2) return "無解";
}
k=(k1+k2)/2;
}
return "查找" + kk + "次未找到解";
}
progText.innerText='擬合函數源程序：\r\n'+nihe;
function test(){
var i, Dx = new Array();
var k =Ck.value-0, b=Cb.value-0, N=Cn.value-0;
for(i=0;i<N;i++) Dx[i] = Math.floor(i*k+b); //序列生成
var k0=Ck0.value-0, b0=Cb0.value-0, a0=Ca0.value-0;
Cjg.value = nihe(Dx, k0, b0, a0);
if(Cjg.value=='ok') Cjg.value = 'k=' + Dx.gk + ' b=' + Dx.gb + ' h=' +Dx.gh;
}
</script>
</td></tr>
</table>
</center>
</body></html>
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -