RCPPARMADILLO FASTLM 결과는 R의 LM과 다르며 무엇을 잘못 했습니까?
문제
다음은 sourceCpp
로 소싱 된 CPP 파일입니다.
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
// [[Rcpp::export]]
List mylm(NumericVector yr, NumericMatrix Xr) {
int n = Xr.nrow(), k = Xr.ncol();
arma::mat X(Xr.begin(), n, k, false); // reuses memory and avoids extra copy
arma::colvec y(yr.begin(), yr.size(), false);
arma::colvec coef = arma::solve(X, y); // fit model y ~ X
arma::colvec resid = y - X*coef; // residuals
double sig2 = arma::as_scalar( arma::trans(resid)*resid/(n-k) );
// std.error of estimate
arma::colvec stderrest = arma::sqrt( sig2 * arma::diagvec( arma::inv(arma::trans(X)*X)) );
return Rcpp::List::create(
Rcpp::Named("coefficients") = coef,
Rcpp::Named("stderr") = stderrest
) ;
}
.
(코드는 @romain françois의 답변 에서이 질문에 이어줍니다 : `sourcecpp 사용``fashlm` 컴파일
그 다음 r :
sourceCpp('b.cpp')
set.seed(1)
x = matrix(rnorm(100), 25, 4)
y = rnorm(25)
mylm(y, x)
summary(lm(y~x))
mylm(y, x)
# $coefficients
# [,1]
# [1,] 0.068978
# [2,] 0.117632
# [3,] -0.029917
# [4,] -0.168648
#
# $stderr
# [,1]
# [1,] 0.17970
# [2,] 0.24312
# [3,] 0.15138
# [4,] 0.21266
summary(lm(y~x))
# Call:
# lm(formula = y ~ x)
#
# Residuals:
# Min 1Q Median 3Q Max
# -0.869 -0.487 -0.271 0.410 1.831
#
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 0.1122 0.1718 0.65 0.52
# x1 0.0499 0.1845 0.27 0.79
# x2 0.1076 0.2470 0.44 0.67
# x3 -0.0435 0.1549 -0.28 0.78
# x4 -0.1750 0.2158 -0.81 0.43
#
# Residual standard error: 0.835 on 20 degrees of freedom
# Multiple R-squared: 0.054, Adjusted R-squared: -0.135
# F-statistic: 0.285 on 4 and 20 DF, p-value: 0.884
. 해결책
기본적으로 R의 LM은 패스를 통과하는 디자인 행렬이 초기 열의 초기 열을 포함하지 않아도 가로 챌 수있는 모델에 적합합니다.따라서 다음은 동일합니다.
lm(y ~ x - 1)
mylm(y, x)
.
'일반'모델을 원하면 디자인 매트릭스를 수정하여 모든 1S의 첫 번째 열을 갖도록해야합니다.
lm(y ~ x)
mylm(y, cbind(1, x))
.
는 동일한 결과를 제공합니다.
제휴하지 않습니다 StackOverflow