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};
线性回归系数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
关于正规方程的注释
首先我不太了解 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};
线性回归系数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
关于正规方程的注释