Merge pull request #1 from donotdespair/main

Week 12 materials from Tomasz

week12/index.qmd

@@ -17,9 +17,27 @@ schedule |>
   pull(ref)
 ```
 
+
 ## What you will learn this week
 
-* Rcpp (with [Tomasz Wozniak](http://bit.ly/tomaszwozniak))
+* **C++** for applications in **R** with [Tomasz Woźniak](https://github.com/donotdespair)
+* The first steps with **Rcpp**
+* Some stats with **RcppArmadillo**
+* Hands-on coding with four exercises
+* Create an **R** package with compiled code in ten steps
+
+## Other references used during the session
+
+* **Armadillo** library [documentation](https://arma.sourceforge.net/docs.html)
+* Tsuda, M., [*Rcpp for everyone*](https://teuder.github.io/rcpp4everyone_en/)
+
+## Scripts used during the session
+
+* Simple list creation [`nicetry.cpp`](https://github.com/numbats/arp/blob/main/week12/nicetry.cpp)
+* Simple linear regression example [`nicelr.cpp`](https://github.com/numbats/arp/blob/main/week12/nicelr.cpp)
+* Simple loop application [`nicelist.cpp`](https://github.com/numbats/arp/blob/main/week12/nicelist.cpp)
+* Sampling from inverted gamma 2 distribution [`nicerig2.cpp`](https://github.com/numbats/arp/blob/main/week12/nicerig2.cpp)
+* Step-by-step package creation [`nicepackage.R`](https://github.com/numbats/arp/blob/main/week12/nicepackage.R)
 
 ```{r}
 #| output: asis

week12/nicelist.cpp

@@ -0,0 +1,12 @@
+#include <Rcpp.h>
+using namespace Rcpp;
+
+// [[Rcpp::export]]
+List nicelist (int n) {
+  NumericVector p = rnorm(n);
+  NumericVector s(n);
+  for (int i=0; i<n; i++) {
+    s[i] =  pow(p[i], 2);
+  }
+  return List::create(_["p"] = p, _["s"] = s);
+}

week12/nicelr.cpp

@@ -0,0 +1,14 @@
+#include <RcppArmadillo.h>
+// [[Rcpp::depends(RcppArmadillo)]]
+using namespace arma;
+
+// [[Rcpp::export]]
+vec nicelr (vec y, mat x) {
+  vec beta_hat = solve(x.t() * x, x.t() * y);
+  return beta_hat;
+}
+
+/*** R
+x = cbind(rep(1,5),1:5); y = x %*% c(1,2) + rnorm(5)
+nicelr(y, x)
+*/

week12/nicepackage.R

@@ -0,0 +1,60 @@
+
+RcppArmadillo::RcppArmadillo.package.skeleton("nicepackage")
+usethis::use_git()
+usethis::use_gpl3_license()
+usethis::use_package_doc()
+roxygen2::roxygenise()
+usethis::use_rcpp_armadillo()
+
+Rcpp::compileAttributes()
+devtools::document()
+devtools::check()
+
+
+devtools::load_all()
+hist(nicerig2(10000, 1, 10), breaks = 100)
+?nicerig2
+?nicepackage
+
+
+
+# to be copy/pasted
+############################################################
+
+#ifndef _NICERIG2_H_
+#define _NICERIG2_H_
+
+#include <RcppArmadillo.h>
+
+arma::vec nicerig2 (
+  const int     n,    // a positive integer - number of draws
+  const double  s,    // a positive scale parameter
+  const double  nu    // a positive shape parameter
+);
+
+#endif  // _NICERIG2_H_
+
+############################################################
+
+#' @title Samples random numbers from the inverted gamma 2 distribution
+#'
+#' @description Provides independent draws
+#' @param n a positive integer, number of draws
+#' @param scale a positive scalar, scale parameter
+#' @param shape a positive integer, shape parameter
+#'
+#' @return a vector with \code{n} independent draws from the inverted gamma 2 distribution
+#'
+#' @export
+
+############################################################
+
+#' @name nicepackage-package
+#' @aliases nicepackage-package nicepackage
+#' @docType package
+#' @importFrom Rcpp sourceCpp
+#' @useDynLib nicepackage, .registration = TRUE
+#' @keywords internal
+"_PACKAGE"
+
+############################################################

week12/nicerig2.cpp

@@ -0,0 +1,17 @@
+#include <RcppArmadillo.h>
+// [[Rcpp::depends(RcppArmadillo)]]
+using namespace arma;
+
+// [[Rcpp::export]]
+vec nicerig2 (
+    const int     n,    // a positive integer - number of draws
+    const double  s,    // a positive scale parameter
+    const double  nu    // a positive shape parameter
+) {
+  vec rig2 = s / chi2rnd( nu, n );
+  return rig2;
+}
+
+/*** R
+nicerig2(2, 1, 1)
+*/

week12/nicetry.cpp

@@ -0,0 +1,15 @@
+#include <Rcpp.h>
+using namespace Rcpp;
+
+// [[Rcpp::export]]
+List nicetry (int n) {
+  NumericVector v = rnorm(n);
+  IntegerVector x = seq_len(n);
+  LogicalVector y = v > 0;
+  CharacterVector z(n, "nice");
+  return List::create(_["v"] = v, _["x"] = x, _["y"] = y, _["z"] = z);
+}
+
+/*** R
+print(nicetry(2))
+*/

week12/nicex1.cpp

@@ -0,0 +1,16 @@
+#include <Rcpp.h>
+using namespace Rcpp;
+
+// [[Rcpp::export]]
+NumericMatrix nicetry (int n) {
+  NumericVector i(n, 1.0);
+  NumericVector t = cumsum(i);
+  NumericVector tt = t - mean(t);
+  NumericVector tt2 = pow(tt, 2);
+  NumericMatrix out = cbind(i, tt, tt2);
+  return out;
+}
+
+/*** R
+nicetry(4)
+*/

week12/nicex2.cpp

@@ -0,0 +1,12 @@
+#include <Rcpp.h>
+using namespace Rcpp;
+
+// [[Rcpp::export]]
+NumericVector nicetry (int n) {
+  NumericVector i = rnorm(n);
+  return cumsum(i);
+}
+
+/*** R
+nicetry(4)
+*/

week12/nicex3.cpp

@@ -0,0 +1,22 @@
+#include <RcppArmadillo.h>
+// [[Rcpp::depends(RcppArmadillo)]]
+
+using namespace Rcpp;
+using namespace arma;
+
+// [[Rcpp::export]]
+List nicelr (vec y, mat x) {
+  mat xx_inv = inv_sympd(x.t() * x);
+  vec beta_hat = xx_inv * x.t() * y;
+  double sigma2 = as_scalar( (y - x * beta_hat).t() * (y - x * beta_hat) / y.n_elem );
+  mat cov_beta_hat = sigma2 * xx_inv;
+  return List::create(
+    _["beta_hat"] = beta_hat,
+    _["cov_beta_hat"] = cov_beta_hat
+  );
+}
+
+/*** R
+x = cbind(rep(1,5),1:5); y = x %*% c(1,2) + rnorm(5)
+nicelr(y, x)
+*/

week12/nicex4.cpp

@@ -0,0 +1,30 @@
+#include <RcppArmadillo.h>
+// [[Rcpp::depends(RcppArmadillo)]]
+
+using namespace Rcpp;
+using namespace arma;
+
+// [[Rcpp::export]]
+List nicernig2 (
+    const int n, 
+    const vec mu,
+    const mat Sigma,
+    const double s, 
+    const double nu
+) {
+
+  vec rig2 = s / chi2rnd( nu, n );
+  mat X(n, mu.n_elem);
+  for (int i=0; i<n; i++) {
+    X.row(i) = trans(mvnrnd( mu, Sigma ));
+  }
+
+  return List::create(
+    _["x"] = X,
+    _["sigma2"] = rig2
+  );
+}
+
+/*** R
+nicernig2(4, rep(0,2), diag(2), 1, 1)
+*/