Python 模块或算法产生与 Mathematica LinearModelFit 相同的结果答案

【问题标题】：Python Module or Algorithm to produce same results as Mathematica LinearModelFitPython 模块或算法产生与 Mathematica LinearModelFit 相同的结果
【发布时间】：2015-09-09 13:24:59
【问题描述】：

首先我并不真正了解 Mathematica，而且我已经很长时间没有做过统计了。

我一直在尝试寻找（Google 和 RTFM）一种方法来使用 scipy.stats.linregress 重现 Mathematica LinearModelFit 函数产生的结果。现在很明显，除了最简单的情况外，这不是要走的路。

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

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

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

提前致谢。

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

Mathematica Notebook 截图：http://imgur.com/owMg3r8

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

【问题讨论】：

您可能想发布实际工作的mathematica 代码（符号名称中不允许使用下划线）和一些示例数据。
我已按要求完成了。
我对mathematica不熟悉，但我建议你看看sklearn.linear_model.LinearRegression。 statsmodels 中也有线性回归功能。我希望这可以帮助

标签： python python-2.7 wolfram-mathematica

【解决方案1】：

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

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

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

相当于：

 xtransform = 1/(2 n - #)^100  & /@ Range[Length[ydata]];
 LinearModelFit[Transpose[{xtransform, ydata}], x, x]["BestFit"]

16.1504 + 1.471945513739138*10^315 x

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

【讨论】：

【解决方案2】：

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

y 数据已导入。

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

线性回归系数a和b是通过求解正规方程得到的。（有关推导，请参见下面的注释）。一旦计算出来，它们就可以在不需要求解器的情况下重复使用。

Clear[a, b, n, Σx, Σy, Σxy, Σx2]

Column[{a, b} = Simplify[First[{a, b} /. Solve[{
    (* Normal equations for straight line *)
    Σy == n a + b Σx,
    Σxy == a Σx + b Σx2},
   {a, b}]]]]

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

【讨论】：