layout: post title: 复杂误差传播的计算 categories:
在数据处理中, 有时我们需要根据拟合参数的误差计算相应变量的误差, 也就是计算误差传播. 对于简单的模型, 计算比较简单, 可以手动完成; 对于复杂的模型, 涉及导数的计算, 处理起来就有点麻烦, 还容易出错. 所以我觉得还是找个比较统一的处理方法才好. 想了一下, 我们需要使用能够进行符号计算的程序, 也就是能够处理公式的程序, 而不仅仅是计算数值. 常用的符号处理程序有mathematica, maple, maxima, matlab, 虽然都能用, 但都太笨重了, 不适合用于处理这么简单的问题. 就想起很老的一个符号计算器EigenMath. 看来一下, 这个程序的Windows版已经不再更新了, 只有一个很老的版本还在网上流传. 不过网站上也提供了源码, 可以自己编译.
下面是一个处理的代码和例子
拟合的模型为
$$y=f(T, P_1, P_2, P_3, P_4)= P_1 - {P_2 \over T} \exp({P_3 \over T^2+P_4})$$
<table class="highlighttable"><th colspan="2" style="text-align:left">error.eigenmath</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height:125%"> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41</pre></div></td><td class="code"><div class="highlight"><pre style="line-height:125%"><span></span><span style="color: #008800"># 定义拟合函数</span> f(T, P1, P2, P3, P4)<span style="color: #666666">=</span> P1 <span style="color: #666666">-</span> (P2<span style="color: #666666">/</span>T)<span style="color: #666666">*</span>exp(P3<span style="color: #666666">/</span>(T<span style="color: #666666">^2+</span>P4)) <span style="color: #008800"># 最终误差的计算根据需要, 可以采用任意一种</span> <span style="color: #008800"># 误差计算方法1: 简单方法</span> f1<span style="color: #666666">=</span>d(f,P1) <span style="#FF0000">#</span> <span style="#FF0000">计算对每个参数的导数</span> f2<span style="color: #666666">=</span>d(f,P2) f3<span style="color: #666666">=</span>d(f,P3) f4<span style="color: #666666">=</span>d(f,P4) df(x,P1,P2,P3,P4 dP1, dP2, dP3, dP4)<span style="color: #666666">=</span>sqrt( (f1<span style="color: #666666">*</span>dP1)<span style="color: #666666">^2+</span>(f2<span style="color: #666666">*</span>dP2)<span style="color: #666666">^2+</span>(f3<span style="color: #666666">*</span>dP3)<span style="color: #666666">^2+</span>(f4<span style="color: #666666">*</span>dP4)<span style="color: #666666">^2</span> ) <span style="color: #008800"># 误差计算方法2: 使用梯度, 数组</span> g<span style="color: #666666">=</span>d(f,(P1, P2, P3, P4)) <span style="#FF0000">#</span> <span style="#FF0000">梯度</span> d<span style="color: #666666">=</span>(dP1, dP2, dP3, dP4) <span style="#FF0000">#</span> <span style="#FF0000">各个参数的误差</span> a<span style="color: #666666">=</span>zero(<span style="color: #666666">4</span>) <span style="color: #AA22FF; font-weight: bold">for</span>(i,<span style="color: #666666">1</span>,<span style="color: #666666">4</span>, a[i]<span style="color: #666666">=</span>g[i]<span style="color: #666666">*</span>d[i]) DeltaF1<span style="color: #666666">=</span>abs(a) DeltaF1 <span style="color: #008800"># 误差计算方法3: 使用数组方法</span> <span style="color: #AA22FF; font-weight: bold">for</span>(i,<span style="color: #666666">1</span>,<span style="color: #666666">4</span>, g[i]<span style="color: #666666">=</span>g[i]<span style="color: #666666">^2</span>, d[i]<span style="color: #666666">=</span>d[i]<span style="color: #666666">^2</span>) <span style="#FF0000">#</span> square DeltaF2<span style="color: #666666">=</span>sqrt(dot(g,d)) DeltaF2 <span style="color: #008800"># 给定参数的拟合值</span> P1<span style="color: #666666">=45.6388</span> P2<span style="color: #666666">=56542</span> P3<span style="color: #666666">=352754</span> P4<span style="color: #666666">=2932.55</span> <span style="color: #008800"># 拟合误差</span> dP1<span style="color: #666666">=0.1013</span> dP2<span style="color: #666666">=65.08</span> dP3<span style="color: #666666">=349</span> dP4<span style="color: #666666">=81.57</span> <span style="color: #008800"># 输出结果</span> <span style="color: #008800">#for(i,0,9,T=400+i*25, print(T,float(df(T,P1,P2,P3,P4 dP1, dP2, dP3, dP4))) )</span> <span style="color: #AA22FF; font-weight: bold">for</span>(i,<span style="color: #666666">0</span>,<span style="color: #666666">1</span>,T<span style="color: #666666">=400+</span>i<span style="color: #666666">*25</span>, print(<span style="color: #00BB00; font-weight: bold">float</span>(DeltaF2)) )</pre></div> </td></tr></table>可以考虑利用js的符号处理库将其改为在线工具, 可以利用的js库比较多了, 随便看到的几个似乎都可以