<p>一个不需要复杂功能的简单算法可以用任何语言编写。在</p>
<p>导入<code>y</code>数据。在</p>
<pre><code>y = {11.56999969, 14.47999954, ... , 340.730011, 202.1699982, 4054.949951};
</code></pre>
<p>线性回归系数<code>a</code>和{<cd3>}通过求解正态方程得到。(推导见下文注释)。一旦计算出来,它们就可以被重用,而不需要求解器。在</p>
^{pr2}$
<blockquote>
<p>(Σx Σxy - Σx2 Σy)/(Σx^2 - n Σx2)</p>
<p>(-n Σxy + Σx Σy)/(Σx^2 - n Σx2)</p>
</blockquote>
<p>根据模型,<code>X</code>被线性化为<code>x</code>。在</p>
<pre><code>n = Length[y]
</code></pre>
<blockquote>
<p>1267</p>
</blockquote>
<pre><code>X = Range[n];
x = Map[1/(2 n - #)^100 &, X];
</code></pre>
<p>计算工程量:</p>
<pre><code>Σx = Sum[x[[i]], {i, n}];
Σy = Sum[y[[i]], {i, n}];
Σxy = Sum[x[[i]]*y[[i]], {i, n}];
Σx2 = Sum[x[[i]]^2, {i, n}];
</code></pre>
<p>执行系数公式:</p>
<pre><code>a = (Σx Σxy - Σx2 Σy)/(Σx^2 - n Σx2)
b = (Σx Σy - n Σxy)/(Σx^2 - n Σx2)
</code></pre>
<blockquote>
<p>16.65767846718208</p>
<p>4.213538401691473*10^313</p>
</blockquote>
<p>在线性化数据上绘制回归线(带缩放)。在</p>
<pre><code>scaled = 10^340;
Show[ListPlot[Transpose[{x scaled, y}],
PlotRange -> {Automatic, {0, 30}}],
ListPlot[Transpose[{x scaled, Table[a + b i, {i, x}]}],
PlotRange -> All, PlotStyle -> Red]]
</code></pre>
<p><a href="https://i.stack.imgur.com/d3jQp.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/d3jQp.png" alt="enter image description here"/></a></p>
<p>重新应用该模型,最小二乘拟合为:<code>a + b/(2 n - X)^100</code></p>
<pre><code>Show[ListPlot[Transpose[{X, y}],
PlotRange -> {Automatic, {0, 400}}],
Plot[a + b/(2 n - X)^100, {X, 0, n},
PlotRange -> {Automatic, {0, 400}}, PlotStyle -> Red]]
</code></pre>
<p><a href="https://i.stack.imgur.com/33IqV.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/33IqV.png" alt="enter image description here"/></a></p>
<p>这与Mathematica的内置解决方案匹配,如下所示。在</p>
<p>同时计算R的平方。在</p>
<pre><code>(* Least-squares regression of y on x *)
Array[(Y[#] = a + b x[[#]]) &, n];
Array[(e[#] = y[[#]] - Y[#]) &, n];
(* Residual or unexplained sum of squares *)
RSS = Sum[e[i]^2, {i, n}];
(* Total sum of squares in the dependent variable, measured about its mean *)
TSS = (y - Mean[y]).(y - Mean[y]);
(* Coefficient of determination, R^2 *)
R2 = 1 - RSS/TSS
</code></pre>
<blockquote>
<p>0.230676</p>
</blockquote>
<p>检查Mathematica的内置功能。在</p>
<pre><code>Clear[x]
lm = LinearModelFit[y, 1/(2 n - x)^100, x];
lm["BestFit"]
</code></pre>
<p><a href="https://i.stack.imgur.com/8oSGW.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/8oSGW.png" alt="enter image description here"/></a></p>
<pre><code>lm["RSquared"]
</code></pre>
<blockquote>
<p>0.230676</p>
</blockquote>
<p><strong>关于正规方程的注记</strong></p>
<p>来源:<a href="https://books.google.com/books/about/Econometric_Methods.html?id=vhKqQgAACAAJ" rel="nofollow noreferrer">Econometric Methods</a></p>
<p><a href="https://i.stack.imgur.com/bsCn6.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/bsCn6.png" alt="enter image description here"/></a></p>