最小二乘法拟合公式的有关分析...

-
--
                               最小二乘法拟合公式的有关分析  
-
                                                                                   史锦顺
-
1 线性拟合公式的应用
       因变量Y与自变量X间，大体有线性关系。
       用实验（实际测量），来确定Y对X的函数关系。测量N次，数据为：
                 X1  X2 …… Xi …… XN  
                 Y1  Y2 …… Yi …… YN                                            （1）
       自变量X取准，误差可略，认为无误差。因变量Y是仪器的测得值，有误差。Yi与Xi一一对应。
       拟合得到的线性公式为：
                  Y = a + bX                                                         （2）
       所谓“最小二乘法拟合”，就是依数据列（1）求得公式（2）。或者说是求得直线方程（2）的斜率b和截距a。
-
       时频界，求频标的漂移率（晶振称老化率）普遍用最小二乘法。
       “自变量”是时间t，ti =iT, T是时间单位（日或半日），误差可略。“因变量”是频标的频率值fi，是仪器测得值，有误差。
       拟合得频率值随时间变化的线性公式为
                  f = fo + Kt                                                          （3）
       将（3）式减去参考频率f标，得频差关系式为：
                  f - f标 = fo - f标 +Kt
                  Δf =Δfo + Kt                                                        （4）
       为方便，通常表示为相对频差
                  Δf /fo =Δfo /fo+ Kt/fo
                  δf = δfo + kt                                                        （5）
-
       要分清两个不同的问题。
       第一个问题是：已知频率变化率k、初始频差偏差δfo 的条件下，怎样用公式（5）求得总偏差δf。
       （5）式是个简单的代数式。δfo与kt都是偏差，它们二者合成，就是代数运算，二者都是有符号、有数值的量，按“代数计算法则”计算就是了，搞不确定度合成，是看错了对象，画蛇添足，瞎胡算。
       例1 已知晶振甲的老化率k= -3×10-10。2017年4月1日，用铯原子钟测得晶振甲的相对频差是 – 1.0×10-8，问到2017年7月10日，长期加电的晶振甲的频差是多少？
       解答：k= -3×10-8/日，δfo= - 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 对晶振甲，已知老化率k= -3×10-10，参照检定规程的频率调整法。2017年4月1日，用铯原子钟调整晶振甲的相对频差是 +1.0×10-8，问到2017年7月10日，长期加电的晶振甲的频差是多少？
       解答：k= -3×10-10/日，δfo= +1.0×10-8，t=100日，代入（5）式，得
                δf = +1.0×10-8 + (-3×10-10/日)×100日
                    = +1.0×10-8  - 3×10-8
                    = - 2×10-8
-
       用最小二乘法，费力气计算得到的公式（5），上两例的应用是正确的。不确定度体系对公式（5）的处理，用什么“不确定度传播律”，是错误的，是抹煞人类的智慧，否定知识。这是“不确定度体系是伪科学”判断的又一个证据。请大家注意一个事实：凡用GUM法处理的计算问题，几乎都是错误的。至于不确定度论说教中的“相关”“不相关”都是没用的假话。初始频差、线性变化率、时间，三者都完全是各自独立的量，该认为是“不相关”；而三者共同构成总的频差的量值，又怎能说不相关？其实初始频差、变化产生的频差，客观上的作用必定是“代数和”，那种关于“相关性”的分析是没有用的。
       如果是估计“范围”，起决定作用的是二项和平方展开式的交叉系数。这里是两项系统误差，算“范围”也必须取“绝对和”，而不能是不确定度体系的“方和根”。“假设不相关”无用；“真不相关”也没用。相关性的判别，与误差合成法无关。
-
       上边讲的是第一个问题：拟合公式的应用；第二个问题是，测定（5）式中两个量截距与斜率的误差有多大？初始频差的测量误差易得；而求老化率k的误差，是当今世界测量计量界的一大难题。叶德培先生在她的样板评定中，缺失了（或者说弄错了）。
       老史先将老化率k的公式简化；由此方见端倪。简述如下。
-
2 老化率公式的简化
       简化公式的主要技巧是对称编号。
       取测量次数N为奇数，中间数为0，则上有(N-1)/2个号，下有(N-1)/2个号。
       例如
       1）测量晶振的频率老化率（线性漂移率），测量7天。每天定时（例如9:00）测量。每日一个数据（三个频率的平均值），对称编号就是f(-3)、f(-2) 、f(-1)、f(0)、f(+1)、f(+2)、f(+3)。
       2）测量晶振的频率老化率。前7天。每天定时两次测量（例如9:00、21:00）。每日两个数据（每个数据是三次测量的平均值），第8天9:00得1数据。对称编号就是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 = fo + Kt                                                          （3）
       其中，fo是初始频率，K是线性变化率。fo易获知；主要作业是拟合直线斜率K。史锦顺得到的简化公式（参见附件）为：
                K = ∑j=(-n)→n (fj –fr) j / [(N-1)N(N+1)T/12]            （6）
       简化公式的主要技巧是对称编号。  
-
【晶振日老化率速算法】
（1） 7天，每天1值
       按时间顺序标记数据：中间数标0，由0分界，顺时序标1至3，逆时序标-1至-3。各数据减一常数之后的尾数乘标号，乘积累加，除以28，即得日老化率。
       口诀：对称编号，去整作差，号乘差累加，除以二十八。  
（2）7周天，每天2点，共15点（第8天测一个点）
       按时间顺序标记数据：中间数标0，由0分界，顺时序标1至7，逆时序标-1至-7。各数据减一常数之后的尾数乘标号，乘积累加，除以140，即得日老化率。
       口诀：对称编号，去整留零，累加号乘零，除以140。
-
3 拟合误差
       拟合误差指求截距与斜率的误差。
       斜率K的简化（严格式）表达为：
                  K = ∑j=(-n)→n (fj –fr)j / [(1/12)(N-1)N(N+1)T]            （6）

       （6）式中，fr可以随意取值，因此K的测量误差与测量fj时的系统误差的恒定值无关。测量晶振老化率用原子频标，其系统误差为恒定值（无频率漂移）。因此误差来自原子频标与比对器的随机误差。随机误差取“方和根”。
       令fr=0,有
                  K =∑j=(-n)→n fj j / [(1/12)(N-1)N(N+1)T]                    （7）
       将fj 表达为：
                  fj = fjo + ξ                                                               （8）
       （8）式代入（7）
                  K =∑j=(-n)→n (fjo+ξ)j / [(1/12)(N-1)N(N+1)T]
                     =∑j=(-n)→n fjo 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）
       测量用同一原子频标与比对器。各次测量，随机误差相同，因此σ可提出来。分子成为j的平方求和，j从1到n，结果乘2.
                  Kσ2=2σ2（12+22+32……+n2）/分母2
                       = 2 (1/6) n(n+1)(2n+1) σ 2 /分母2 （查数学手册得知自然数平方之和）
                       = (1/12)(N-1)N(N+1) σ 2  / 分母2
-
                  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）
-
       原子频标的随机误差σ到老化率的传递系数，Kσ /σ
       取样间隔时间：日
               N=3    1/√2
               N=5    1/√10
               N=7    1/√28
       取样间隔时间：0.5日（12小时）
               N=15   1/√70
-
-
附件 论最小二乘法
-
论最小二乘法.doc
-
包括有关最小二乘法计算的三个材料：
《新概念测量学》（中卷）p8
《新概念测量学》（下卷）p52
《驳不确定度论一百六十篇集》p340
-

你说的很对。就像史老师说的，如果是整个模型来看的话，当然相关了。
但是我还是关注于y=ax+b这个直线的基础模型。a和x都是输入量么？应该说是的，而且a确实也是由x得出的。
但是a在这里并非是由x直接得出的，而是根据最小二乘法得出的。得出的a是斜率，一个定值。一旦计算出来，a是不变的，而x是一个随机变量。
所以并非说是一个简单的勾股定理的问题，或者说用啥勾股定理也不对，这里是斜率，就是用也是用a=y/x，前天我脑子短路了。不过，我觉得直接用这个也不对，最小二乘法中间还有其他步骤，而不是直接得出斜率。

如果我写：
y=ax+b
a=y/x
上面的方程组来看，因为y在第一个式子作伪输出量，第二个式子作伪输入量，应该是属于一个负反馈模型，别的我不知道，我学过负反馈的作用是产生一个恒稳态，也就是一个定值，难道就是说这个a么？我觉得两个式子放在一起就是一个非线性的模型了。
就像史老师说的k值，虽然t会变化，但是因为分子的f把其约掉了，所以最后的k是不变的，两者还是不相关。



这里我觉得以下的几个结论不知道对不对：

如果一个输入量是另外一个输入量得出的，相关的可能性很大。但是不见得必然相关，比如上面这个情况，要具体问题具体分析才行。

-
1 分子分母的对消作用
       你只注意到分母中有T（采样周期），忽视了分子中频率变化。一个周期（T），频率变化是KT, 则2个周期是2KT，三个周期是3KT，分子与分母中T的作用对消了，因此测得值K与T无关。
       T取1日时K的测得值，同T=0.5日时K的测得值，若测量设备（原子钟与比对器）误差可略时，二者是相同的。由于频差测量的误差，T值取大，则一个周期的频差大，易于测准。但又要考虑工作量。对新研制的原子频标（如铷频标）的漂移率的测量，T取一星期，则需42周日（N=7）。T取1个月则需整半年（6周月）。晶振是通用仪器，不可能太费时。
-
2 相关性的对象
       不确定度体系（包括某些现代误差理论书籍）所讲究的“相关性”，指的是参与合成的各个量的误差的相关性。不是函数与自变量之间的相关性。函数与自变量之间有明确的函数关系，不存在“不相关”的问题。
       上例之对消情况是特例。求得的K值与T在函数的层次上不相关，但通常函数与自变量是相关的。误差合成（不确定度合成）的所谓相关性，是指分母的T与分子的频率偏差值是否相关。很明显二者强相关。这是在“频差值”这个层面上说的。但在频差值的测量误差（原子频标与比对器引入）这个层面上说，有另有相关性问题。频率值的误差与时间的测定误差，又显然不相关。
-
3 讲究相关性，是误导
       相关系数的概念来自统计理论，是针对随机变量而言的，对测量计量中的随机误差，可用。但对系统误差，相关系数的概念没用。
       1）相关系数的公式，仅能用于随机误差。对二系统误差，计算相关系数仅能得零。就是说，既有随机误差又有系统误差的场合，没有能用于判断相关性的公式。
       2） 不确定度评定，都“假设不相关”，那是掩耳盗铃。
       3） GUM与JJF上有关“有系统误差就可忽略协方差”的条款都是误导，是错误的。
-
4 客观起作用的是交叉系数
       误差合成（不确定度合成）法，取决于交叉系数。由此，不管相关不相关，两项大系统误差合成必须取“绝对和”。
-
5 不确定度体系的合成之路走不通
       B类不确定度评定，认定仪器误差范围是“均匀分布”，统计方法错位。对时域统计，系统误差是恒值，是窄脉冲分布。系统误差无方差可言。
       二系统误差项的相关系数不能求。假设“不相关”，是反科学的。科学技术不能假设。
-
       总之当今的不确定度合成，对系统误差，分布搞错，相关系数不能求，是条走不通的死路。
-

附件的内容才是真详实，爽。说实话内容对我这种不怎么懂时间测量的，要啃有些吃力，有的看了。
不过我不太明白，想请教史老师，线性变化率，也就是∑j=(-n)→n (fj –fr) j / [(N-1)N(N+1)T/12]  这个，T为采样周期，说明线性变化率是采样周期算出来的，T的改变会引起K的改变。可是为什么说线性变化率和时间是独立的量呢？

y=a+bx 我知道您在考虑之前的那个问题，但我感觉那个公式中sa，sb，sx进行合成，是因为这个a和b在用最小二乘法计算时，其实求出的是一个最佳拟合线，故a和b存在标准差，或者说存在一个范围的。就好比不确定度评定中，一个标准电阻阻值R，它必然是个定值，但它会依然引入一个不确定度分量。
R=U/I，按这个公式R是由U和I计算出来的，但大家都知道电阻阻值R是定值，是和U，I无关的。。。如果真要计算相关性的话，应该相关性为0.

非常感谢史老师的解读。若是没你的提醒，我还觉得频率等于时间的倒数，看来我错把这里的周期等同于时间了。
前段时间我还在问一个关于最小二乘法的问题，现在我觉得可以从几何的更一般意义上去说：
y=a+bx这个式子，b是斜率，是一个定值，x是一个变化的量，x怎样变化b都不改变，所以他们本身应该是独立的
不知道这样理解可对？

对，我突然想到，x怎么变，a不是还是那个值吗？前段时间我考虑用勾股定理这种想法得出相关的结论，是不是有待商榷
其实就是真的一个值是另外一个值得出的话，他们一定相关么？

对。理解正确。

对。就是这样。

其实就是真的一个值是另外一个值得出的话，他们一定相关么？

个人感觉,一个值是另外一个值得出的话，那一个是输入量，一个是输出量，它们之间有对应的函数关系，就不应该用相关不相关来表述了。按照不确定度评定中来看。两输入量之间的关系用相关性，相关系数表述。而输入量和输出量之间得关系用函数式和灵敏系数表示。