生成与Mathematica LinearModelFi相同结果的Python模块或算法

Question

首先，我对Mathematica不太了解，而且我已经很长时间没有做过统计了。在

我一直在尝试（Google和RTFM）使用scipy.stats.linregress来重现MathematicaLinearModelFit函数产生的结果。现在很明显，除了最简单的情况外，这是不可能的。在

LinearModelFit[ydata, 1/(2 n - x)^100, x]

产生16.3766 + <<70>>/(2580 - x)^100

如果有人能为我指出正确的方向，我将不胜感激。在

提前谢谢。在

数据：http://pastebin.com/RTp5em0W

Mathematica笔记本的屏幕截图：http://imgur.com/owMg3r8

注：我没有做数学工作。ddd是可以在pastebin链接上找到的数据。分母中的y应该是x

Answer 1

一个不需要复杂功能的简单算法可以用任何语言编写。在

导入y数据。在

y = {11.56999969, 14.47999954, ... , 340.730011, 202.1699982, 4054.949951};

线性回归系数a和{}通过求解正态方程得到。（推导见下文注释）。一旦计算出来，它们就可以被重用，而不需要求解器。在

^{pr2}$

(Σx Σxy - Σx2 Σy)/(Σx^2 - n Σx2)

(-n Σxy + Σx Σy)/(Σx^2 - n Σx2)

根据模型，X被线性化为x。在

n = Length[y]

1267

X = Range[n];
x = Map[1/(2 n - #)^100 &, X];

计算工程量：

Σ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}];

执行系数公式：

a = (Σx Σxy - Σx2 Σy)/(Σx^2 - n Σx2)
b = (Σx Σy - n Σxy)/(Σx^2 - n Σx2)

16.65767846718208

4.213538401691473*10^313

在线性化数据上绘制回归线（带缩放）。在

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]]

重新应用该模型，最小二乘拟合为：a + b/(2 n - X)^100

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]]

这与Mathematica的内置解决方案匹配，如下所示。在

同时计算R的平方。在

(* 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

0.230676

检查Mathematica的内置功能。在

Clear[x]

lm = LinearModelFit[y, 1/(2 n - x)^100, x];
lm["BestFit"]

lm["RSquared"]

0.230676

关于正规方程的注记

来源：Econometric Methods

Answer 2

我不知道python的解决方案，但解决这个问题的一种方法是根据作为LinearModelFit参数提供的函数形式转换x数据：

 n=1290
 LinearModelFit[ydata, 1/(2 n - x)^100, x]["BestFit"]

16.1504 + 1.471945513739138*10^315/(2580 - x)^100

相当于：

^{pr2}$

16.1504 + 1.471945513739138*10^315 x

您应该能够很容易地进行这种转换，并在python中使用标准线性回归。但是，由于指数较大，您可能会遇到精度问题。在