数学不好的人就是郁闷啊。
前两天写的自动外理ＨＰＬＣ数据的程序当中有一段代码是用来计算标准曲线的斜率的。这东西简单啊，一元线性回归嘛，用最小二乘法，以前写过的。
void Min2Method(float &slope, float &intercept, double X[], double Y[], int nCount) { double SumX, SumY, SumXY, SumX2; SumX = 0; SumX2 = 0; for(int i=0; i<nCount; i++) { SumX += X[i]; SumX2 += (X[i] * X[i]); } SumY = 0; for(int i=0; i<nCount; i++) { SumY += Y[i]; } SumXY = 0; for(int i=0; i<nCount; i++) { SumXY += ( X[i] * Y[i]); } intercept = ( (SumX2*SumY - SumX*SumXY) / (nCount*SumX2 - SumX * SumX)); slope = ( (nCount*SumXY - SumX*SumY) / (nCount*SumX2 - SumX * SumX)); }
可是，老板要求，标准曲线要过原点，也就是过零点的斜率。要把横轴截距设置为零。这下就又不会了。只好google research（因为手头没有书）。谁让咱统计学没学好呢。
经过一番学习，终于了重新找回了那么一点点感觉。首先还是搞清楚上面的代码是怎么来的吧。
设有一组输入变量：（x1, x2, … , xn）都对应的有一组输出变量：（y1,y2, … , yn），它们满足一次线性关系：
y = a + b * x
其中 yi = a + b * xi + εi ( i = 1, …, n); 这里的期望值E(ε)=0, 方差Var(ε)＝σ2属随机误差。
最小二乘法是高斯发明的，找到一条直线, 确定a和b的值使偏差平方和 (Q = Σ[yi - (a+b * xi)]2) 最小，即它对a和b求偏导为零：
(∂Q/∂a)= -2Σ[yi - (a+b * xi)] = 0 (1)
(∂Q/∂b)= -2Σ{xi * [yi - (a+b * xi)]} = 0 (2)
分别得解：
n * a + b * Σ(xi) = Σ(yi) (3)
a * Σ(xi)+ b * Σ(xi2) = Σ(xiyi) (4)
得到a和b
b = Σ[(xi-Σ(xi)/n) * (yi-Σ(yi)/n)]/Σ[(xi-Σ(xi)/n)2];
a = [Σ(yi) - b * Σ(xi)] / n.
现在如果要求a为零，那就是说（1）式不成立，只需要求解(2)式就可以了。当a=0, (2)变化为：
(∂Q/∂b)= -2Σ[xi * (yi - b * xi)] = 0 （5）
解得：
b * Σ(xi2) = Σ(xiyi)
b = Σ(xiyi)/Σ(xi2)
于是上面的程序为了得到过原点的直线斜率就可以改写为：
void Min2Method(float &slope, float &R2, double X[], double Y[], int nCount) { double SumX, SumY, SumXY, SumX2, SumY2, SumSSE, SumSST; SumX = 0; SumX2 = 0; for(int i=0; i<nCount; i++) { SumX += X[i]; SumX2 += (X[i] * X[i]); } SumY = 0; SumY2 = 0; for(int i=0; i<nCount; i++) { SumY += Y[i]; SumY2 += (Y[i] * Y[i]); } SumXY = 0; for(int i=0; i<nCount; i++) { SumXY += ( X[i] * Y[i]); } slope = SumXY / SumX2; SumSSE = 0; double avgY = SumY/nCount; for(int i=0; i<nCount; i++) { SumSSE += (Y[i] - avgY)*(Y[i] - avgY); SumSST += (Y[i] - X[i]*slope) * (Y[i] - X[i]*slope); } R2 = 1 - SumSST/SumSSE; }
转载请注明文章来自糗世界博客
呵呵，我手头还有高数的书，上面有最小二乘法。不过现在都直接用Excel算了。
[回复]
admin 回复:
三月 4th, 2010 at 7:13 下午
@Liang Huang, 用Excel就是慢啊……如果一个星期要处理20条标准曲线，然后还要从上百个样品原文件里读取面积而后算浓度，两三个小时别的不干，就光copy paste了。
[回复]