最小二乘法擬合公式的有關(guān)分析

-
--
                               最小二乘法擬合公式的有關(guān)分析  
-
                                                                                   史錦順
-
1 線性擬合公式的應(yīng)用
       因變量Y與自變量X間，大體有線性關(guān)系。
       用實(shí)驗(yàn)（實(shí)際測量），來確定Y對(duì)X的函數(shù)關(guān)系。測量N次，數(shù)據(jù)為：
                 X₁  X₂ …… X_i …… X_N
                 Y₁  Y₂ …… Y_i …… Y_N                                            （1）
       自變量X取準(zhǔn)，誤差可略，認(rèn)為無誤差。因變量Y是儀器的測得值，有誤差。Y_i與X_i一一對(duì)應(yīng)。
       擬合得到的線性公式為：
                  Y = a + bX                                                         （2）
       所謂“最小二乘法擬合”，就是依數(shù)據(jù)列（1）求得公式（2）。或者說是求得直線方程（2）的斜率b和截距a。
-
       時(shí)頻界，求頻標(biāo)的漂移率（晶振稱老化率）普遍用最小二乘法。
       “自變量”是時(shí)間t，t_i =iT, T是時(shí)間單位（日或半日），誤差可略。“因變量”是頻標(biāo)的頻率值fi，是儀器測得值，有誤差。
       擬合得頻率值隨時(shí)間變化的線性公式為
                  f = f_o + Kt                                                          （3）
       將（3）式減去參考頻率f_標(biāo)，得頻差關(guān)系式為：
                  f - f_標(biāo) = f_o - f_標(biāo) +Kt
                  Δf =Δf_o + Kt                                                        （4）
       為方便，通常表示為相對(duì)頻差
                  Δf /f_o =Δf_o /f_o+ Kt/f_o
                  δf = δf_o + kt                                                        （5）
-
       要分清兩個(gè)不同的問題。
       第一個(gè)問題是：已知頻率變化率k、初始頻差偏差δfo 的條件下，怎樣用公式（5）求得總偏差δf。
       （5）式是個(gè)簡單的代數(shù)式。δf_o與kt都是偏差，它們二者合成，就是代數(shù)運(yùn)算，二者都是有符號(hào)、有數(shù)值的量，按“代數(shù)計(jì)算法則”計(jì)算就是了，搞不確定度合成，是看錯(cuò)了對(duì)象，畫蛇添足，瞎胡算。
       例1 已知晶振甲的老化率k= -3×10^-10。2017年4月1日，用銫原子鐘測得晶振甲的相對(duì)頻差是 – 1.0×10^-8，問到2017年7月10日，長期加電的晶振甲的頻差是多少？
       解答：k= -3×10^-8/日，δf_o= - 1.0×10^-8，t=100日，代入（5）式，得
                δf = - 1.0×10^-8 + (-3×10^-10/日)×100日
                    = - 1.0×10^-8 - 3×10^-8
                    = - 4×10^-8
-
       例2 對(duì)晶振甲，已知老化率k= -3×10^-10，參照檢定規(guī)程的頻率調(diào)整法。2017年4月1日，用銫原子鐘調(diào)整晶振甲的相對(duì)頻差是 +1.0×10^-8，問到2017年7月10日，長期加電的晶振甲的頻差是多少？
       解答：k= -3×10-¹⁰/日，δf_o= +1.0×10^-8，t=100日，代入（5）式，得
                δf = +1.0×10^-8 + (-3×10-¹⁰/日)×100日
                    = +1.0×10^-8  - 3×10^-8
                    = - 2×10^-8
-
       用最小二乘法，費(fèi)力氣計(jì)算得到的公式（5），上兩例的應(yīng)用是正確的。不確定度體系對(duì)公式（5）的處理，用什么“不確定度傳播律”，是錯(cuò)誤的，是抹煞人類的智慧，否定知識(shí)。這是“不確定度體系是偽科學(xué)”判斷的又一個(gè)證據(jù)。請(qǐng)大家注意一個(gè)事實(shí)：凡用GUM法處理的計(jì)算問題，幾乎都是錯(cuò)誤的。至于不確定度論說教中的“相關(guān)”“不相關(guān)”都是沒用的假話。初始頻差、線性變化率、時(shí)間，三者都完全是各自獨(dú)立的量，該認(rèn)為是“不相關(guān)”；而三者共同構(gòu)成總的頻差的量值，又怎能說不相關(guān)？其實(shí)初始頻差、變化產(chǎn)生的頻差，客觀上的作用必定是“代數(shù)和”，那種關(guān)于“相關(guān)性”的分析是沒有用的。
       如果是估計(jì)“范圍”，起決定作用的是二項(xiàng)和平方展開式的交叉系數(shù)。這里是兩項(xiàng)系統(tǒng)誤差，算“范圍”也必須取“絕對(duì)和”，而不能是不確定度體系的“方和根”。“假設(shè)不相關(guān)”無用；“真不相關(guān)”也沒用。相關(guān)性的判別，與誤差合成法無關(guān)。
-
       上邊講的是第一個(gè)問題：擬合公式的應(yīng)用；第二個(gè)問題是，測定（5）式中兩個(gè)量截距與斜率的誤差有多大？初始頻差的測量誤差易得；而求老化率k的誤差，是當(dāng)今世界測量計(jì)量界的一大難題。葉德培先生在她的樣板評(píng)定中，缺失了（或者說弄錯(cuò)了）。
       老史先將老化率k的公式簡化；由此方見端倪。簡述如下。
-
2 老化率公式的簡化
       簡化公式的主要技巧是對(duì)稱編號(hào)。
       取測量次數(shù)N為奇數(shù)，中間數(shù)為0，則上有(N-1)/2個(gè)號(hào)，下有(N-1)/2個(gè)號(hào)。
       例如
       1）測量晶振的頻率老化率（線性漂移率），測量7天。每天定時(shí)（例如9:00）測量。每日一個(gè)數(shù)據(jù)（三個(gè)頻率的平均值），對(duì)稱編號(hào)就是f(-3)、f(-2) 、f(-1)、f(0)、f(+1)、f(+2)、f(+3)。
       2）測量晶振的頻率老化率。前7天。每天定時(shí)兩次測量（例如9:00、21:00）。每日兩個(gè)數(shù)據(jù)（每個(gè)數(shù)據(jù)是三次測量的平均值），第8天9:00得1數(shù)據(jù)。對(duì)稱編號(hào)就是f(-7)、f(-6) 、f(-5)、f(-4)、f(-3)、f(-2)、f(-1)、f(0)、f(+1) 、f(+2)、f(+3)、f(+4)、f(+5) 、f(+6)、f(+7)。
-
       要擬合的公式是
                 f = f_o + Kt                                                          （3）
       其中，fo是初始頻率，K是線性變化率。fo易獲知；主要作業(yè)是擬合直線斜率K。史錦順得到的簡化公式（參見附件）為：
                K = ∑_j=(-n)→n (f_j –f_r) j / [(N-1)N(N+1)T/12]            （6）
       簡化公式的主要技巧是對(duì)稱編號(hào)。  
-
【晶振日老化率速算法】
（1） 7天，每天1值
       按時(shí)間順序標(biāo)記數(shù)據(jù)：中間數(shù)標(biāo)0，由0分界，順時(shí)序標(biāo)1至3，逆時(shí)序標(biāo)-1至-3。各數(shù)據(jù)減一常數(shù)之后的尾數(shù)乘標(biāo)號(hào)，乘積累加，除以28，即得日老化率。
       口訣：對(duì)稱編號(hào)，去整作差，號(hào)乘差累加，除以二十八。  
（2）7周天，每天2點(diǎn)，共15點(diǎn)（第8天測一個(gè)點(diǎn)）
       按時(shí)間順序標(biāo)記數(shù)據(jù)：中間數(shù)標(biāo)0，由0分界，順時(shí)序標(biāo)1至7，逆時(shí)序標(biāo)-1至-7。各數(shù)據(jù)減一常數(shù)之后的尾數(shù)乘標(biāo)號(hào)，乘積累加，除以140，即得日老化率。
       口訣：對(duì)稱編號(hào)，去整留零，累加號(hào)乘零，除以140。
-
3 擬合誤差
       擬合誤差指求截距與斜率的誤差。
       斜率K的簡化（嚴(yán)格式）表達(dá)為：
                  K = ∑_j=(-n)→n (f_j –f_r)j / [(1/12)(N-1)N(N+1)T]            （6）

       （6）式中，f_r可以隨意取值，因此K的測量誤差與測量fj時(shí)的系統(tǒng)誤差的恒定值無關(guān)。測量晶振老化率用原子頻標(biāo)，其系統(tǒng)誤差為恒定值（無頻率漂移）。因此誤差來自原子頻標(biāo)與比對(duì)器的隨機(jī)誤差。隨機(jī)誤差取“方和根”。
       令f_r=0,有
                  K =∑_j=(-n)→n f_j j / [(1/12)(N-1)N(N+1)T]                    （7）
       將f_j 表達(dá)為：
                  f_j = f_jo + ξ                                                               （8）
       （8）式代入（7）
                  K =∑_j=(-n)→n (f_jo+ξ)j / [(1/12)(N-1)N(N+1)T]
                     =∑_j=(-n)→n f_jo j / [(1/12)(N-1)N(N+1)T]
                        +∑_j=(-n)→n ξ j / [(1/12)(N-1)N(N+1)T]                  （9）
       由（9）知，K 的誤差為
                  ΔK=∑_j=(-n)→n ξ j / [(1/12)(N-1)N(N+1)T]
                  K_σ=∑_j=(-n)→n σ j / [(1/12)(N-1)N(N+1)T]                     （10）
       測量用同一原子頻標(biāo)與比對(duì)器。各次測量，隨機(jī)誤差相同，因此σ可提出來。分子成為j的平方求和，j從1到n，結(jié)果乘2.
                  K_σ²=2σ²（1²+2²+3²……+n²）/分母2
                       = 2 (1/6) n(n+1)(2n+1) σ ² /分母² （查數(shù)學(xué)手冊得知自然數(shù)平方之和）
                       = (1/12)(N-1)N(N+1) σ ² / 分母²
-
                  K_σ = σ√[(1/12)(N-1)N(N+1)] / 分母
                      = σ√[(1/12)(N-1)N(N+1)] / [(1/12)(N-1)N(N+1)T]

                  K_σ =(σ/ T) /√[(1/12)(N-1)N(N+1)]                            （11）
-
       原子頻標(biāo)的隨機(jī)誤差σ到老化率的傳遞系數(shù)，K_σ /σ
       取樣間隔時(shí)間：日
               N=3    1/√2
               N=5    1/√10
               N=7    1/√28
       取樣間隔時(shí)間：0.5日（12小時(shí)）
               N=15   1/√70
-
-
附件 論最小二乘法
-
論最小二乘法.doc (292 KB, 下載次數(shù): 46)

2017-4-21 11:53 上傳

點(diǎn)擊文件名下載附件
下載積分: 金幣 -1

-
包括有關(guān)最小二乘法計(jì)算的三個(gè)材料：
《新概念測量學(xué)》（中卷）p8
《新概念測量學(xué)》（下卷）p52
《駁不確定度論一百六十篇集》p340
-

附件的內(nèi)容才是真詳實(shí)，爽。說實(shí)話內(nèi)容對(duì)我這種不怎么懂時(shí)間測量的，要啃有些吃力，有的看了。
不過我不太明白，想請(qǐng)教史老師，線性變化率，也就是∑j=(-n)→n (fj –fr) j / [(N-1)N(N+1)T/12]  這個(gè)，T為采樣周期，說明線性變化率是采樣周期算出來的，T的改變會(huì)引起K的改變。可是為什么說線性變化率和時(shí)間是獨(dú)立的量呢？

solarup 發(fā)表于 2017-4-21 13:04
附件的內(nèi)容才是真詳實(shí)，爽。說實(shí)話內(nèi)容對(duì)我這種不怎么懂時(shí)間測量的，要啃有些吃力，有的看了。
不過我不太 ...

-
1 分子分母的對(duì)消作用
       你只注意到分母中有T（采樣周期），忽視了分子中頻率變化。一個(gè)周期（T），頻率變化是KT, 則2個(gè)周期是2KT，三個(gè)周期是3KT，分子與分母中T的作用對(duì)消了，因此測得值K與T無關(guān)。
       T取1日時(shí)K的測得值，同T=0.5日時(shí)K的測得值，若測量設(shè)備（原子鐘與比對(duì)器）誤差可略時(shí)，二者是相同的。由于頻差測量的誤差，T值取大，則一個(gè)周期的頻差大，易于測準(zhǔn)。但又要考慮工作量。對(duì)新研制的原子頻標(biāo)（如銣頻標(biāo)）的漂移率的測量，T取一星期，則需42周日（N=7）。T取1個(gè)月則需整半年（6周月）。晶振是通用儀器，不可能太費(fèi)時(shí)。
-
2 相關(guān)性的對(duì)象
       不確定度體系（包括某些現(xiàn)代誤差理論書籍）所講究的“相關(guān)性”，指的是參與合成的各個(gè)量的誤差的相關(guān)性。不是函數(shù)與自變量之間的相關(guān)性。函數(shù)與自變量之間有明確的函數(shù)關(guān)系，不存在“不相關(guān)”的問題。
       上例之對(duì)消情況是特例。求得的K值與T在函數(shù)的層次上不相關(guān)，但通常函數(shù)與自變量是相關(guān)的。誤差合成（不確定度合成）的所謂相關(guān)性，是指分母的T與分子的頻率偏差值是否相關(guān)。很明顯二者強(qiáng)相關(guān)。這是在“頻差值”這個(gè)層面上說的。但在頻差值的測量誤差（原子頻標(biāo)與比對(duì)器引入）這個(gè)層面上說，有另有相關(guān)性問題。頻率值的誤差與時(shí)間的測定誤差，又顯然不相關(guān)。
-
3 講究相關(guān)性，是誤導(dǎo)
       相關(guān)系數(shù)的概念來自統(tǒng)計(jì)理論，是針對(duì)隨機(jī)變量而言的，對(duì)測量計(jì)量中的隨機(jī)誤差，可用。但對(duì)系統(tǒng)誤差，相關(guān)系數(shù)的概念沒用。
       1）相關(guān)系數(shù)的公式，僅能用于隨機(jī)誤差。對(duì)二系統(tǒng)誤差，計(jì)算相關(guān)系數(shù)僅能得零。就是說，既有隨機(jī)誤差又有系統(tǒng)誤差的場合，沒有能用于判斷相關(guān)性的公式。
       2） 不確定度評(píng)定，都“假設(shè)不相關(guān)”，那是掩耳盜鈴。
       3） GUM與JJF上有關(guān)“有系統(tǒng)誤差就可忽略協(xié)方差”的條款都是誤導(dǎo)，是錯(cuò)誤的。
-
4 客觀起作用的是交叉系數(shù)
       誤差合成（不確定度合成）法，取決于交叉系數(shù)。由此，不管相關(guān)不相關(guān)，兩項(xiàng)大系統(tǒng)誤差合成必須取“絕對(duì)和”。
-
5 不確定度體系的合成之路走不通
       B類不確定度評(píng)定，認(rèn)定儀器誤差范圍是“均勻分布”，統(tǒng)計(jì)方法錯(cuò)位。對(duì)時(shí)域統(tǒng)計(jì)，系統(tǒng)誤差是恒值，是窄脈沖分布。系統(tǒng)誤差無方差可言。
       二系統(tǒng)誤差項(xiàng)的相關(guān)系數(shù)不能求。假設(shè)“不相關(guān)”，是反科學(xué)的。科學(xué)技術(shù)不能假設(shè)。
-
       總之當(dāng)今的不確定度合成，對(duì)系統(tǒng)誤差，分布搞錯(cuò)，相關(guān)系數(shù)不能求，是條走不通的死路。
-

史錦順 發(fā)表于 2017-4-22 12:30
-
1 分子分母的對(duì)消作用
       你只注意到分母中有T（采樣周期），忽視了分子中頻率變化。一個(gè)周期（T ...

非常感謝史老師的解讀。若是沒你的提醒，我還覺得頻率等于時(shí)間的倒數(shù)，看來我錯(cuò)把這里的周期等同于時(shí)間了。
前段時(shí)間我還在問一個(gè)關(guān)于最小二乘法的問題，現(xiàn)在我覺得可以從幾何的更一般意義上去說：
y=a+bx這個(gè)式子，b是斜率，是一個(gè)定值，x是一個(gè)變化的量，x怎樣變化b都不改變，所以他們本身應(yīng)該是獨(dú)立的
不知道這樣理解可對(duì)？

solarup 發(fā)表于 2017-4-24 16:36
非常感謝史老師的解讀。若是沒你的提醒，我還覺得頻率等于時(shí)間的倒數(shù)，看來我錯(cuò)把這里的周期等同于時(shí)間了 ...

是不是意思的， X1  X2 …… Xi …… XN,Y1  Y2 …… Yi …… YN   給定后，a和b根據(jù)最小二乘法求解，就是個(gè)定值了，y=a+bx，而之后x無論如何的選取，都對(duì)這個(gè)函數(shù)式無任何的影響了？

吳下阿蒙 發(fā)表于 2017-4-24 17:15
是不是意思的， X1  X2 …… Xi …… XN,Y1  Y2 …… Yi …… YN   給定后，a和b根據(jù)最小二乘法求解，就 ...

對(duì)，我突然想到，x怎么變，a不是還是那個(gè)值嗎？前段時(shí)間我考慮用勾股定理這種想法得出相關(guān)的結(jié)論，是不是有待商榷
其實(shí)就是真的一個(gè)值是另外一個(gè)值得出的話，他們一定相關(guān)么？

solarup 發(fā)表于 2017-4-24 16:36
非常感謝史老師的解讀。若是沒你的提醒，我還覺得頻率等于時(shí)間的倒數(shù)，看來我錯(cuò)把這里的周期等同于時(shí)間了 ...

       對(duì)。理解正確。

吳下阿蒙 發(fā)表于 2017-4-24 17:15
是不是意思的， X1  X2 …… Xi …… XN,Y1  Y2 …… Yi …… YN   給定后，a和b根據(jù)最小二乘法求解，就 ...

       對(duì)。就是這樣。

solarup 發(fā)表于 2017-4-25 07:38
對(duì)，我突然想到，x怎么變，a不是還是那個(gè)值嗎？前段時(shí)間我考慮用勾股定理這種想法得出相關(guān)的結(jié)論，是不是 ...

其實(shí)就是真的一個(gè)值是另外一個(gè)值得出的話，他們一定相關(guān)么？

個(gè)人感覺,一個(gè)值是另外一個(gè)值得出的話，那一個(gè)是輸入量，一個(gè)是輸出量，它們之間有對(duì)應(yīng)的函數(shù)關(guān)系，就不應(yīng)該用相關(guān)不相關(guān)來表述了。按照不確定度評(píng)定中來看。兩輸入量之間的關(guān)系用相關(guān)性，相關(guān)系數(shù)表述。而輸入量和輸出量之間得關(guān)系用函數(shù)式和靈敏系數(shù)表示。

吳下阿蒙 發(fā)表于 2017-4-25 09:17
其實(shí)就是真的一個(gè)值是另外一個(gè)值得出的話，他們一定相關(guān)么？

個(gè)人感覺,一個(gè)值是另外一個(gè)值得出的話，那 ...

你說的很對(duì)。就像史老師說的，如果是整個(gè)模型來看的話，當(dāng)然相關(guān)了。
但是我還是關(guān)注于y=ax+b這個(gè)直線的基礎(chǔ)模型。a和x都是輸入量么？應(yīng)該說是的，而且a確實(shí)也是由x得出的。
但是a在這里并非是由x直接得出的，而是根據(jù)最小二乘法得出的。得出的a是斜率，一個(gè)定值。一旦計(jì)算出來，a是不變的，而x是一個(gè)隨機(jī)變量。
所以并非說是一個(gè)簡單的勾股定理的問題，或者說用啥勾股定理也不對(duì)，這里是斜率，就是用也是用a=y/x，前天我腦子短路了。不過，我覺得直接用這個(gè)也不對(duì)，最小二乘法中間還有其他步驟，而不是直接得出斜率。

如果我寫：
y=ax+b
a=y/x
上面的方程組來看，因?yàn)閥在第一個(gè)式子作偽輸出量，第二個(gè)式子作偽輸入量，應(yīng)該是屬于一個(gè)負(fù)反饋模型，別的我不知道，我學(xué)過負(fù)反饋的作用是產(chǎn)生一個(gè)恒穩(wěn)態(tài)，也就是一個(gè)定值，難道就是說這個(gè)a么？我覺得兩個(gè)式子放在一起就是一個(gè)非線性的模型了。
就像史老師說的k值，雖然t會(huì)變化，但是因?yàn)榉肿拥膄把其約掉了，所以最后的k是不變的，兩者還是不相關(guān)。



這里我覺得以下的幾個(gè)結(jié)論不知道對(duì)不對(duì)：

如果一個(gè)輸入量是另外一個(gè)輸入量得出的，相關(guān)的可能性很大。但是不見得必然相關(guān)，比如上面這個(gè)情況，要具體問題具體分析才行。

solarup 發(fā)表于 2017-4-25 10:37
你說的很對(duì)。就像史老師說的，如果是整個(gè)模型來看的話，當(dāng)然相關(guān)了。
但是我還是關(guān)注于y=ax+b這個(gè)直線的 ...

y=a+bx 我知道您在考慮之前的那個(gè)問題，但我感覺那個(gè)公式中sa，sb，sx進(jìn)行合成，是因?yàn)檫@個(gè)a和b在用最小二乘法計(jì)算時(shí)，其實(shí)求出的是一個(gè)最佳擬合線，故a和b存在標(biāo)準(zhǔn)差，或者說存在一個(gè)范圍的。就好比不確定度評(píng)定中，一個(gè)標(biāo)準(zhǔn)電阻阻值R，它必然是個(gè)定值，但它會(huì)依然引入一個(gè)不確定度分量。
R=U/I，按這個(gè)公式R是由U和I計(jì)算出來的，但大家都知道電阻阻值R是定值，是和U，I無關(guān)的。。。如果真要計(jì)算相關(guān)性的話，應(yīng)該相關(guān)性為0.

吳下阿蒙 發(fā)表于 2017-4-25 12:09
y=a+bx 我知道您在考慮之前的那個(gè)問題，但我感覺那個(gè)公式中sa，sb，sx進(jìn)行合成，是因?yàn)檫@個(gè)a和b在用最小 ...

我覺得，要考慮模型和現(xiàn)實(shí)之間的差距
歐姆定律算出來的R是定值，和最小二乘法的斜率是定值，兩者是不同的。
最小二乘法中的斜率定值，是x是隨機(jī)變量，但是無論在概率上取何值，擬合出來的直線斜率是不變的。
但是歐姆定律不是，U和I不同時(shí)，R的值并非定值，U和I在概率上取值的不同，會(huì)得出不同的R，只有他們并非隨機(jī)變量時(shí)，才能得出固定的R值

非常感謝。

非常感謝。

吳下阿蒙 發(fā)表于 2017-4-25 12:09
y=a+bx 我知道您在考慮之前的那個(gè)問題，但我感覺那個(gè)公式中sa，sb，sx進(jìn)行合成，是因?yàn)檫@個(gè)a和b在用最小 ...

是有個(gè)"應(yīng)用場景"問題要事先明確---

如果"系統(tǒng)參數(shù)"a、b由某一組"輸入"、"輸出"觀測數(shù)據(jù)[x]~[y]"擬合"所得，那它們與此組[x]的"相關(guān)性"通常應(yīng)該是"不容忽視"的。……但實(shí)際應(yīng)用中可能不存在通過琢磨其"相關(guān)性"的大小來"估算"那組"輸出"數(shù)據(jù)[y]的"散布"大小---這組[y]的觀測數(shù)據(jù)是"已知"的，直接"統(tǒng)計(jì)"可得其"散布"的大小，弄清爽了不會(huì)"繞"著轉(zhuǎn)圈。

若是"擬合"所用以外的另一組"輸入"數(shù)據(jù)[x]，則擬合所得的"系統(tǒng)"參數(shù)與其"相關(guān)性"通常或許是可以忽略的。

這個(gè)必須要留言，保留了學(xué)習(xí)，謝謝計(jì)量前輩無私的分享~~！！！！