想改进这个问题吗?Update the question所以堆栈溢出的值小于aa>。
两年前关闭。
我试图将线性回归R结果与python的结果相匹配
匹配每个自变量的系数
下面是代码:
数据已上载。
on-topic
https://www.dropbox.com/s/oowe4irm9332s78/X.csv?dl=0
R代码:
#define pathname = " "
X <- read.csv(file.path(pathname,"X.csv"),stringsAsFactors = F)
Y <- read.csv(file.path(pathname,"Y.csv"),stringsAsFactors = F)
d1 = Y$d1
reg_data <- cbind(d1, X)
head(reg_data)
reg_model <- lm(d1 ~ .,data=reg_data)
names(which(is.na(reg_model$coefficients)))
reg_model$coefficients
R结果
> summary(reg_model)
Call:
lm(formula = d1 ~ ., data = reg_data)
Residuals:
ALL 60 residuals are 0: no residual degrees of freedom!
Coefficients: (14 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) -67.37752 NA NA NA
v1 1.30214 NA NA NA
v2 -2.93118 NA NA NA
v3 7.58902 NA NA NA
v4 11.88570 NA NA NA
v5 1.60622 NA NA NA
v6 3.71528 NA NA NA
v7 -9.34627 NA NA NA
v8 -3.84694 NA NA NA
v9 -2.51332 NA NA NA
v10 4.22403 NA NA NA
v11 -9.70126 NA NA NA
v12 NA NA NA NA
v13 4.67276 NA NA NA
v14 -6.57924 NA NA NA
v15 -3.68065 NA NA NA
v16 5.25168 NA NA NA
v17 14.60444 NA NA NA
v18 16.00679 NA NA NA
v19 24.79622 NA NA NA
v20 13.85774 NA NA NA
v21 2.16022 NA NA NA
v22 -36.65361 NA NA NA
v23 2.26554 NA NA NA
v24 NA NA NA NA
v25 NA NA NA NA
v26 7.00981 NA NA NA
v27 0.88904 NA NA NA
v28 0.34400 NA NA NA
v29 -5.27597 NA NA NA
v30 5.21034 NA NA NA
v31 6.79640 NA NA NA
v32 2.96346 NA NA NA
v33 -1.52702 NA NA NA
v34 -2.74632 NA NA NA
v35 -2.36952 NA NA NA
v36 -7.76547 NA NA NA
v37 2.19630 NA NA NA
v38 1.63336 NA NA NA
v39 0.69485 NA NA NA
v40 0.37379 NA NA NA
v41 -0.09107 NA NA NA
v42 2.06569 NA NA NA
v43 1.57505 NA NA NA
v44 2.70535 NA NA NA
v45 1.17634 NA NA NA
v46 -10.51141 NA NA NA
v47 -1.15060 NA NA NA
v48 2.87353 NA NA NA
v49 3.37740 NA NA NA
v50 -5.89816 NA NA NA
v51 0.85851 NA NA NA
v52 3.73929 NA NA NA
v53 4.93265 NA NA NA
v54 3.45650 NA NA NA
v55 0.12382 NA NA NA
v56 -0.21171 NA NA NA
v57 4.37199 NA NA NA
v58 3.21456 NA NA NA
v59 0.09012 NA NA NA
v60 -0.85414 NA NA NA
v61 -3.29856 NA NA NA
v62 4.38842 NA NA NA
v63 NA NA NA NA
v64 NA NA NA NA
v65 NA NA NA NA
v66 NA NA NA NA
v67 NA NA NA NA
v68 NA NA NA NA
v69 NA NA NA NA
v70 NA NA NA NA
v71 NA NA NA NA
v72 NA NA NA NA
v73 NA NA NA NA
Residual standard error: NaN on 0 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: NaN
F-statistic: NaN on 59 and 0 DF, p-value: NA
Python代码:
Y = pd.read_csv(pathname+"Y.csv")
X = pd.read_csv(pathname+"X.csv")
lr = LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)
lr.fit(X, Y['d1'])
(list(zip(lr.coef_, X)))
lr.intercept_
Python结果:
intercept = 29.396033164254106
[(-2.4463986167806304, 'v1'),
(-1.6293010275307021, 'v2'),
(0.89089949009506508, 'v3'),
(-3.1021251646895251, 'v4'),
(-1.7707078771936109, 'v5'),
(-2.0474705122225636, 'v6'),
(-1.5537181337496202, 'v7'),
(-1.6391241229716156, 'v8'),
(-1.2981646048517046, 'v9'),
(0.89221826294889328, 'v10'),
(-0.56694104645951571, 'v11'),
(2.042810365310288e-14, 'v12'),
(-2.0312478672439052, 'v13'),
(-1.5617121392788413, 'v14'),
(0.4583365939498274, 'v15'),
(0.8840538748922967, 'v16'),
(-5.5952681002058871, 'v17'),
(2.4937042448512892, 'v18'),
(0.45806845189176543, 'v19'),
(-1.1648810657830406, 'v20'),
(-1.7800004329275585, 'v21'),
(-5.0132817522704816, 'v22'),
(3.6862778096189266, 'v23'),
(2.7533531010703882e-14, 'v24'),
(1.2150003225741557e-14, 'v25'),
(0.94669823515018103, 'v26'),
(-0.3082823207975679, 'v27'),
(0.53619247380957358, 'v28'),
(-1.1339902793546781, 'v29'),
(1.9657159583080186, 'v30'),
(-0.63200501460653324, 'v31'),
(1.4741013580918978, 'v32'),
(-2.4448418291953313, 'v33'),
(-2.0787115960875036, 'v34'),
(0.22492914212063603, 'v35'),
(-0.75136276693004922, 'v36'),
(1.2838658951186761, 'v37'),
(0.5816277993227944, 'v38'),
(-0.11270569554555088, 'v39'),
(-0.13430982360936233, 'v40'),
(-3.3189296496897662, 'v41'),
(-0.452575588270415, 'v42'),
(6.1329755709937519, 'v43'),
(0.61559185634634817, 'v44'),
(-1.206315459828555, 'v45'),
(-3.7452010299772009, 'v46'),
(-1.1472174665136678, 'v47'),
(2.8960489381172652, 'v48'),
(0.0090220136972478659, 'v49'),
(-5.264918363314754, 'v50'),
(1.2194758337662015, 'v51'),
(2.78655271320092, 'v52'),
(3.106513852668896, 'v53'),
(3.5181252502607929, 'v54'),
(-0.34426523770507278, 'v55'),
(-0.48792823932479878, 'v56'),
(0.12284460490031779, 'v57'),
(1.6860388628044991, 'v58'),
(1.2823067194737174, 'v59'),
(2.8352263554153665, 'v60'),
(-1.304336378501032, 'v61'),
(0.55226132316435139, 'v62'),
(1.5416988124754771, 'v63'),
(-0.2605804175310813, 'v64'),
(1.2489066081702334, 'v65'),
(-0.44469553013696161, 'v66'),
(-1.4102990055550157, 'v67'),
(3.8150423259022639, 'v68'),
(0.12039684410168072, 'v69'),
(-1.340699466779357, 'v70'),
(1.7066389124439212, 'v71'),
(0.50470752944860442, 'v72'),
(1.0024872633969766, 'v73')]
但不匹配。请帮忙。
注:以下示例与之匹配
https://www.dropbox.com/s/79scp54unzlbwyk/Y.csv?dl=0
最佳答案
dr如果你想在Python中复制R的策略,你可能需要自己实现它,因为R做了一些在其他地方不常见的聪明的事情。
作为参考(目前仅在评论中提到),这是一种不适定/秩亏拟合,当预测变量多于响应时(p>n:在这种情况下p=73,n=61),通常当存在许多分类响应和/或实验设计在某种程度上受到限制时,会发生这种情况。处理这些情况以便得到一个答案,这意味着任何事情通常都需要仔细思考和/或先进的技术(例如,惩罚回归:参见Wikipedia article on linear regression中的套索和脊回归参考文献)。
处理这种情况最幼稚的方法是将其全部投入到标准的线性代数中,并希望不会出现太严重的中断,这显然就是python的statsmodels包所做的:从pitfalls document开始:
秩亏矩阵不会产生误差。
几乎完全多重共线性或病态设计矩阵的情况可能会产生数值不稳定的结果。如果这不是所需的行为,用户需要手动检查矩阵的秩或条件数。
下一个最好的事情(当共线度很小时合理)是在做线性代数时明智地旋转,也就是说,重新安排计算问题,这样共线部分就可以省略掉。这就是R所做的;为了做到这一点,R代码的作者必须修改标准的LINPACK例程。
underlying code有以下评论/解释:
c dqrdc2使用householder变换计算qr
用p矩阵x.a有限列对n的c因式分解
基于缩减列2-范数的c轴策略
c将具有接近零范数的列移到
c x矩阵。这个策略意味着
自由度效应可以用自然的方法计算。
我对用这种方式修改linpack代码感到非常紧张。
如果你是计算线性代数大师
c了解如何解决这个问题请随时
c建议对此代码进行改进。
这段代码(和注释)从1998年起就在R的代码库中;我很想知道最初是谁编写的(基于comments further down in the code它似乎是Ross Ihaka?),但除了1998年的代码重组之外,我很难跟踪代码的历史。(稍微深入的研究表明,这段代码基本上是自R的记录历史开始以来就在R的代码库中,也就是说,该文件是在SVN第2版1997-09-18中添加的,直到后来才被修改。)
Martin Mächler最近(2016年10月25日,here)添加了更多关于?qr
的信息,因此这些信息将在下一版本的R。。。
如果你知道如何将编译好的FORTRAN代码与Python代码链接起来(我不知道),那么编译src/appl/dqrdc2.f
并将lm.fit
的内核转换成Python将非常容易:这是lm.fit
的核心,减去错误检查和其他处理。。。
z <- .Call(C_Cdqrls, x, y, tol, FALSE)
coef <- z$coefficients
pivot <- z$pivot
r1 <- seq_len(z$rank)
dn <- colnames(x)
nmeffects <- c(dn[pivot[r1]], rep.int("", n - z$rank))
r2 <- if (z$rank < p)
(z$rank + 1L):p
else integer()
if (is.matrix(y)) { ## ...
} else {
coef[r2] <- NA
if (z$pivoted)
coef[pivot] <- coef
names(coef) <- dn
names(z$effects) <- nmeffects
}
z$coefficients <- coef
关于python - R和Python中线性回归的差异,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/40935624/