View difference between Paste ID: 2xJbCjdd and eMM8yuWF

 library(kernlab)
 N = 10000
-sourceCpp(code = '#include <Rcpp.h>
+p = 5
 b = 1000
 X = matrix(rnorm(N*p), ncol = p)
 center = X[sample(1:N, b),]
-DataFrame mypnorm(NumericVector x){
+sigma = 1
-	int n = x.size();
+kernel_X <- kernelMatrix(rbfdot(sigma=1/(2*sigma^2)), X, center)
-	NumericVector y1(n), y2(n), y3(n);
+
-
+library(Rcpp)
-	for (int i=0; i<n; i++){
+library(RcppArmadillo)
-		y1[i] = ::Rf_pnorm5(x[i], 0.0, 1.0, 1, 0);
+
-		y2[i] = R::pnorm(x[i], 0.0, 1.0, 1, 0);
+
 // [[Rcpp::depends(RcppArmadillo)]]
-	y3 = pnorm(x);
+
-	return DataFrame::create(Named("R") = y1,
+#include <omp.h>
-				 Named("Rf_") = y2,
+
-				 Named("sugar") = y3);
+using namespace arma;
 // [[Rcpp::export]]
-mypnorm(runif(10, -3, 3))
+
 NumericMatrix kernelMatrix_cpp(NumericMatrix Xr, NumericMatrix Centerr, double sigma) {
 	omp_set_num_threads(omp_get_max_threads());
 	uword n = Xr.nrow(), b = Centerr.nrow(), row_index, col_index;
 	mat X(Xr.begin(), n, Xr.ncol(), false);
 	mat Center(Centerr.begin(), b, Centerr.ncol(), false);
 	mat KerX(n, b);
 	#pragma omp parallel private(row_index, col_index)
-List fastLm(NumericVector yr, NumericMatrix Xr) {
+	for (row_index = 0; row_index < n; row_index++)
-	int n = Xr.nrow(), k = Xr.ncol();
+	{
-	arma::mat X(Xr.begin(), n, k, false);
+		#pragma omp for nowait
-	arma::colvec y(yr.begin(), yr.size(), false);
+		for (col_index = 0; col_index < b; col_index++)
-	arma::colvec coef = arma::solve(X, y);
+		{
-	arma::colvec resid = y - X*coef;
+			KerX(row_index, col_index) = exp(sum(square(X.row(row_index) - Center.row(col_index))) / (-2.0 * sigma * sigma));
-	double sig2 = arma::as_scalar(arma::trans(resid)*resid/(n-k));
+		}
-	arma::colvec stderrest = arma::sqrt(sig2 * arma::diagvec( arma::inv(arma::trans(X)*X)) );
+
-	return List::create(Named("coefficients") = coef,
+    return wrap(KerX);
-                            Named("stderr") = stderrest);
+
 kernel_X_cpp = kernelMatrix_cpp(X, center, sigma)
-X = matrix(rnorm(100*2), ncol = 2)
+all.equal([email protected], kernel_X_cpp)
-y = 0.5 + X %*% c(0.3, -0.7) + rnorm(100)
+sum(abs([email protected] - kernel_X_cpp) < 1e-12)
-fastLm(y, X)
+
 library(rbenchmark)
-sourceCpp(code = '// [[Rcpp::depends(RcppEigen)]]
+benchmark(kernelMatrix_cpp(X, center, sigma), kernelMatrix(rbfdot(sigma=1/(2*sigma^2)), X, center), columns=c("test", "replications","elapsed", "relative"), replications=50, order="relative")