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

Python Module or Algorithm to produce same results as Mathematica LinearModelFit

首先我不太了解 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 笔记本的屏幕截图:http://imgur.com/owMg3r8

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

我不知道 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 中进行转换并使用标准线性回归。但是,由于指数很大,您可能会遇到精度问题。

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

已导入 y 数据。

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

线性回归系数ab是通过求解正规方程得到的。 (有关推导,请参见下面的注释)。一旦计算出来,它们就可以重复使用而无需求解器。

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