RcppNumerical: Numerical integration and optimization with Rcpp

[This article was first published on Yixuan's Blog - R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Introduction

I have seen several conversations in Rcpp-devel mailing list asking how to
compute numerical integration or optimization in Rcpp. While R in fact
has the functions Rdqags, Rdqagi, nmmin, vmmin etc. in its API
to accomplish such tasks, it is not so straightforward to use them with Rcpp.

For my own research projects I need to do a lot of numerical integration,
root finding and optimization, so to make my life a little bit easier, I
just created the RcppNumerical
package that simplifies these procedures. I haven’t submitted RcppNumerical
to CRAN, since the API may change quickly according to my needs or the
feedbacks from other people.

Basically RcppNumerical includes a number of open source libraries for
numerical computing, so that Rcpp code can link to this package to use
the functions provided by these libraries. Alternatively, RcppNumerical
provides some wrapper functions that have less configuration and fewer
arguments, if you just want to use the default and quickly get the results.

RcppNumerical depends on Rcpp (obviously) and RcppEigen,

  • To use RcppNumerical with Rcpp::sourceCpp(), add the following two lines
    to the C++ source file:
<span class="c1">// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
</span>
  • To use RcppNumerical in your package, add the corresponding fields to the
    DESCRIPTION file:
Imports: RcppNumerical
LinkingTo: Rcpp, RcppEigen, RcppNumerical

Numerical Integration

(Picture from Wikipedia)

The numerical integration code contained in RcppNumerical is based
on the NumericalIntegration
library developed by Sreekumar Thaithara Balan,
Mark Sauder, and Matt Beall.

To compute integration of a function, first define a functor inherited from
the Func class:

<span class="k">class</span> <span class="nc">Func</span>
<span class="p">{</span>
<span class="k">public</span><span class="o">:</span>
    <span class="k">virtual</span> <span class="kt">double</span> <span class="k">operator</span><span class="p">()(</span><span class="k">const</span> <span class="kt">double</span><span class="o">&</span> <span class="n">x</span><span class="p">)</span> <span class="k">const</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span>
    <span class="k">virtual</span> <span class="kt">void</span>   <span class="k">operator</span><span class="p">()(</span><span class="kt">double</span><span class="o">*</span> <span class="n">x</span><span class="p">,</span> <span class="k">const</span> <span class="kt">int</span> <span class="n">n</span><span class="p">)</span> <span class="k">const</span>
    <span class="p">{</span>
        <span class="k">for</span><span class="p">(</span><span class="kt">int</span> <span class="n">i</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span> <span class="n">i</span> <span class="o"><</span> <span class="n">n</span><span class="p">;</span> <span class="n">i</span><span class="o">++</span><span class="p">)</span>
            <span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="k">this</span><span class="o">-></span><span class="k">operator</span><span class="p">()(</span><span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">]);</span>
    <span class="p">}</span>
<span class="p">};</span>

The first function evaluates one point at a time, and the second version
overwrites each point in the array by the corresponding function values.
Only the second function will be used by the integration code, but usually it
is easier to implement the first one.

RcppNumerical provides a wrapper function for the NumericalIntegration
library with the following interface:

<span class="kr">inline</span> <span class="kt">double</span> <span class="n">integrate</span><span class="p">(</span>
    <span class="k">const</span> <span class="n">Func</span><span class="o">&</span> <span class="n">f</span><span class="p">,</span> <span class="k">const</span> <span class="kt">double</span><span class="o">&</span> <span class="n">lower</span><span class="p">,</span> <span class="k">const</span> <span class="kt">double</span><span class="o">&</span> <span class="n">upper</span><span class="p">,</span>
    <span class="kt">double</span><span class="o">&</span> <span class="n">err_est</span><span class="p">,</span> <span class="kt">int</span><span class="o">&</span> <span class="n">err_code</span><span class="p">,</span>
    <span class="k">const</span> <span class="kt">int</span> <span class="n">subdiv</span> <span class="o">=</span> <span class="mi">100</span><span class="p">,</span>
    <span class="k">const</span> <span class="kt">double</span><span class="o">&</span> <span class="n">eps_abs</span> <span class="o">=</span> <span class="mf">1e-8</span><span class="p">,</span> <span class="k">const</span> <span class="kt">double</span><span class="o">&</span> <span class="n">eps_rel</span> <span class="o">=</span> <span class="mf">1e-6</span><span class="p">,</span>
    <span class="k">const</span> <span class="n">Integrator</span><span class="o"><</span><span class="kt">double</span><span class="o">>::</span><span class="n">QuadratureRule</span> <span class="n">rule</span> <span class="o">=</span> <span class="n">Integrator</span><span class="o"><</span><span class="kt">double</span><span class="o">>::</span><span class="n">GaussKronrod41</span>
<span class="p">)</span>

See the README page for the
explanation of each argument. Below shows an example that calculates the
moment generating function of a $Beta(a,b)$ distribution,
$M(t) = E(e^{tX})$:

<span class="c1">// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
</span><span class="cp">#include <RcppNumerical.h>
</span><span class="k">using</span> <span class="k">namespace</span> <span class="n">Numer</span><span class="p">;</span>

<span class="c1">// M(t) = E(exp(t * X)) = int exp(t * x) * f(x) dx, f(x) is the p.d.f.
</span><span class="k">class</span> <span class="nc">Mintegrand</span><span class="o">:</span> <span class="k">public</span> <span class="n">Func</span>
<span class="p">{</span>
<span class="k">private</span><span class="o">:</span>
    <span class="k">const</span> <span class="kt">double</span> <span class="n">a</span><span class="p">;</span>
    <span class="k">const</span> <span class="kt">double</span> <span class="n">b</span><span class="p">;</span>
    <span class="k">const</span> <span class="kt">double</span> <span class="n">t</span><span class="p">;</span>
<span class="k">public</span><span class="o">:</span>
    <span class="n">Mintegrand</span><span class="p">(</span><span class="kt">double</span> <span class="n">a_</span><span class="p">,</span> <span class="kt">double</span> <span class="n">b_</span><span class="p">,</span> <span class="kt">double</span> <span class="n">t_</span><span class="p">)</span> <span class="o">:</span> <span class="n">a</span><span class="p">(</span><span class="n">a_</span><span class="p">),</span> <span class="n">b</span><span class="p">(</span><span class="n">b_</span><span class="p">),</span> <span class="n">t</span><span class="p">(</span><span class="n">t_</span><span class="p">)</span> <span class="p">{}</span>

    <span class="kt">double</span> <span class="k">operator</span><span class="p">()(</span><span class="k">const</span> <span class="kt">double</span><span class="o">&</span> <span class="n">x</span><span class="p">)</span> <span class="k">const</span>
    <span class="p">{</span>
        <span class="k">return</span> <span class="n">std</span><span class="o">::</span><span class="n">exp</span><span class="p">(</span><span class="n">t</span> <span class="o">*</span> <span class="n">x</span><span class="p">)</span> <span class="o">*</span> <span class="n">R</span><span class="o">::</span><span class="n">dbeta</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">a</span><span class="p">,</span> <span class="n">b</span><span class="p">,</span> <span class="mi">0</span><span class="p">);</span>
    <span class="p">}</span>
<span class="p">};</span>

<span class="c1">// [[Rcpp::export]]
</span><span class="kt">double</span> <span class="nf">beta_mgf</span><span class="p">(</span><span class="kt">double</span> <span class="n">t</span><span class="p">,</span> <span class="kt">double</span> <span class="n">a</span><span class="p">,</span> <span class="kt">double</span> <span class="n">b</span><span class="p">)</span>
<span class="p">{</span>
    <span class="n">Mintegrand</span> <span class="n">f</span><span class="p">(</span><span class="n">a</span><span class="p">,</span> <span class="n">b</span><span class="p">,</span> <span class="n">t</span><span class="p">);</span>
    <span class="kt">double</span> <span class="n">err_est</span><span class="p">;</span>
    <span class="kt">int</span> <span class="n">err_code</span><span class="p">;</span>
    <span class="k">return</span> <span class="n">integrate</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">,</span> <span class="mf">1.0</span><span class="p">,</span> <span class="n">err_est</span><span class="p">,</span> <span class="n">err_code</span><span class="p">);</span>
<span class="p">}</span>

We can compile and run this code in R and draw the graph:

<span class="n">library</span><span class="p">(</span><span class="n">Rcpp</span><span class="p">)</span>
<span class="n">library</span><span class="p">(</span><span class="n">ggplot2</span><span class="p">)</span>
<span class="n">sourceCpp</span><span class="p">(</span><span class="s2">"somefile.cpp"</span><span class="p">)</span>
<span class="n">t0</span> <span class="o">=</span> <span class="n">seq</span><span class="p">(</span><span class="m">-3</span><span class="p">,</span> <span class="m">3</span><span class="p">,</span> <span class="n">by</span> <span class="o">=</span> <span class="m">0.1</span><span class="p">)</span>
<span class="n">mt</span> <span class="o">=</span> <span class="n">sapply</span><span class="p">(</span><span class="n">t0</span><span class="p">,</span> <span class="n">beta_mgf</span><span class="p">,</span> <span class="n">a</span> <span class="o">=</span> <span class="m">1</span><span class="p">,</span> <span class="n">b</span> <span class="o">=</span> <span class="m">1</span><span class="p">)</span>
<span class="n">qplot</span><span class="p">(</span><span class="n">t0</span><span class="p">,</span> <span class="n">mt</span><span class="p">,</span> <span class="n">geom</span> <span class="o">=</span> <span class="s2">"line"</span><span class="p">,</span> <span class="n">xlab</span> <span class="o">=</span> <span class="s2">"t"</span><span class="p">,</span> <span class="n">ylab</span> <span class="o">=</span> <span class="s2">"M(t)"</span><span class="p">,</span>
      <span class="n">main</span> <span class="o">=</span> <span class="s2">"Moment generating function of Beta(1, 1)"</span><span class="p">)</span>

Numerical Optimization

Currently RcppNumerical contains the L-BFGS algorithm for unconstrained
minimization problems based on the
libLBFGS library
developed by Naoaki Okazaki.

(Picture from aria42.com)

Again, one needs to first define a functor to represent the multivariate
function to be minimized.

<span class="k">class</span> <span class="nc">MFuncGrad</span>
<span class="p">{</span>
<span class="k">public</span><span class="o">:</span>
    <span class="k">virtual</span> <span class="kt">double</span> <span class="n">f_grad</span><span class="p">(</span><span class="n">Constvec</span><span class="o">&</span> <span class="n">x</span><span class="p">,</span> <span class="n">Refvec</span> <span class="n">grad</span><span class="p">)</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span>
<span class="p">};</span>

Here Constvec represents a read-only vector and Refvec a writable
vector. Their definitions are

<span class="c1">// Reference to a vector
</span><span class="k">typedef</span> <span class="n">Eigen</span><span class="o">::</span><span class="n">Ref</span><span class="o"><</span><span class="n">Eigen</span><span class="o">::</span><span class="n">VectorXd</span><span class="o">></span>             <span class="n">Refvec</span><span class="p">;</span>
<span class="k">typedef</span> <span class="k">const</span> <span class="n">Eigen</span><span class="o">::</span><span class="n">Ref</span><span class="o"><</span><span class="k">const</span> <span class="n">Eigen</span><span class="o">::</span><span class="n">VectorXd</span><span class="o">></span> <span class="n">Constvec</span><span class="p">;</span>

(Basically you can treat Refvec as a Eigen::VectorXd and
Constvec the const version. Using Eigen::Ref is mainly to avoid
memory copy. See the explanation
here.)

The f_grad() member function returns the function value on vector x,
and overwrites grad by the gradient.

The wrapper function for libLBFGS is

<span class="kr">inline</span> <span class="kt">int</span> <span class="n">optim_lbfgs</span><span class="p">(</span>
    <span class="n">MFuncGrad</span><span class="o">&</span> <span class="n">f</span><span class="p">,</span> <span class="n">Refvec</span> <span class="n">x</span><span class="p">,</span> <span class="kt">double</span><span class="o">&</span> <span class="n">fx_opt</span><span class="p">,</span>
    <span class="k">const</span> <span class="kt">int</span> <span class="n">maxit</span> <span class="o">=</span> <span class="mi">300</span><span class="p">,</span>
    <span class="k">const</span> <span class="kt">double</span><span class="o">&</span> <span class="n">eps_f</span> <span class="o">=</span> <span class="mf">1e-6</span><span class="p">,</span> <span class="k">const</span> <span class="kt">double</span><span class="o">&</span> <span class="n">eps_g</span> <span class="o">=</span> <span class="mf">1e-5</span>
<span class="p">)</span>

Also refer to the README page for
details and see the logistic regression example below.

Fast Logistic Regression: An Example

Let’s see a realistic example that uses the optimization library to fit a
logistic regression.

Given a data matrix $X$ and a 0-1 valued vector $Y$, we want to find a
coefficient vector $beta$ such that the negative log-likelihood function is
minimized:

The gradient function is

So we can write the code as follows:

<span class="c1">// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
</span><span class="cp">#include <RcppNumerical.h>
</span><span class="k">using</span> <span class="k">namespace</span> <span class="n">Numer</span><span class="p">;</span>
<span class="k">using</span> <span class="n">Rcpp</span><span class="o">::</span><span class="n">NumericVector</span><span class="p">;</span>
<span class="k">using</span> <span class="n">Rcpp</span><span class="o">::</span><span class="n">NumericMatrix</span><span class="p">;</span>
<span class="k">typedef</span> <span class="n">Eigen</span><span class="o">::</span><span class="n">Map</span><span class="o"><</span><span class="n">Eigen</span><span class="o">::</span><span class="n">MatrixXd</span><span class="o">></span> <span class="n">MapMat</span><span class="p">;</span>
<span class="k">typedef</span> <span class="n">Eigen</span><span class="o">::</span><span class="n">Map</span><span class="o"><</span><span class="n">Eigen</span><span class="o">::</span><span class="n">VectorXd</span><span class="o">></span> <span class="n">MapVec</span><span class="p">;</span>

<span class="k">class</span> <span class="nc">LogisticReg</span><span class="o">:</span> <span class="k">public</span> <span class="n">MFuncGrad</span>
<span class="p">{</span>
<span class="k">private</span><span class="o">:</span>
    <span class="k">const</span> <span class="n">MapMat</span> <span class="n">X</span><span class="p">;</span>
    <span class="k">const</span> <span class="n">MapVec</span> <span class="n">Y</span><span class="p">;</span>
<span class="k">public</span><span class="o">:</span>
    <span class="n">LogisticReg</span><span class="p">(</span><span class="k">const</span> <span class="n">MapMat</span> <span class="n">x_</span><span class="p">,</span> <span class="k">const</span> <span class="n">MapVec</span> <span class="n">y_</span><span class="p">)</span> <span class="o">:</span> <span class="n">X</span><span class="p">(</span><span class="n">x_</span><span class="p">),</span> <span class="n">Y</span><span class="p">(</span><span class="n">y_</span><span class="p">)</span> <span class="p">{}</span>

    <span class="kt">double</span> <span class="n">f_grad</span><span class="p">(</span><span class="n">Constvec</span><span class="o">&</span> <span class="n">beta</span><span class="p">,</span> <span class="n">Refvec</span> <span class="n">grad</span><span class="p">)</span>
    <span class="p">{</span>
        <span class="c1">// Negative log likelihood
</span>        <span class="c1">//   sum(log(1 + exp(X * beta))) - y' * X * beta
</span>
        <span class="n">Eigen</span><span class="o">::</span><span class="n">VectorXd</span> <span class="n">xbeta</span> <span class="o">=</span> <span class="n">X</span> <span class="o">*</span> <span class="n">beta</span><span class="p">;</span>
        <span class="k">const</span> <span class="kt">double</span> <span class="n">yxbeta</span> <span class="o">=</span> <span class="n">Y</span><span class="p">.</span><span class="n">dot</span><span class="p">(</span><span class="n">xbeta</span><span class="p">);</span>
        <span class="c1">// X * beta => exp(X * beta)
</span>        <span class="n">xbeta</span> <span class="o">=</span> <span class="n">xbeta</span><span class="p">.</span><span class="n">array</span><span class="p">().</span><span class="n">exp</span><span class="p">();</span>
        <span class="k">const</span> <span class="kt">double</span> <span class="n">f</span> <span class="o">=</span> <span class="p">(</span><span class="n">xbeta</span><span class="p">.</span><span class="n">array</span><span class="p">()</span> <span class="o">+</span> <span class="mf">1.0</span><span class="p">).</span><span class="n">log</span><span class="p">().</span><span class="n">sum</span><span class="p">()</span> <span class="o">-</span> <span class="n">yxbeta</span><span class="p">;</span>

        <span class="c1">// Gradient
</span>        <span class="c1">//   X' * (p - y), p = exp(X * beta) / (1 + exp(X * beta))
</span>
        <span class="c1">// exp(X * beta) => p
</span>        <span class="n">xbeta</span><span class="p">.</span><span class="n">array</span><span class="p">()</span> <span class="o">/=</span> <span class="p">(</span><span class="n">xbeta</span><span class="p">.</span><span class="n">array</span><span class="p">()</span> <span class="o">+</span> <span class="mf">1.0</span><span class="p">);</span>
        <span class="n">grad</span><span class="p">.</span><span class="n">noalias</span><span class="p">()</span> <span class="o">=</span> <span class="n">X</span><span class="p">.</span><span class="n">transpose</span><span class="p">()</span> <span class="o">*</span> <span class="p">(</span><span class="n">xbeta</span> <span class="o">-</span> <span class="n">Y</span><span class="p">);</span>

        <span class="k">return</span> <span class="n">f</span><span class="p">;</span>
    <span class="p">}</span>
<span class="p">};</span>

<span class="c1">// [[Rcpp::export]]
</span><span class="n">NumericVector</span> <span class="nf">logistic_reg</span><span class="p">(</span><span class="n">NumericMatrix</span> <span class="n">x</span><span class="p">,</span> <span class="n">NumericVector</span> <span class="n">y</span><span class="p">)</span>
<span class="p">{</span>
    <span class="k">const</span> <span class="n">MapMat</span> <span class="n">xx</span> <span class="o">=</span> <span class="n">Rcpp</span><span class="o">::</span><span class="n">as</span><span class="o"><</span><span class="n">MapMat</span><span class="o">></span><span class="p">(</span><span class="n">x</span><span class="p">);</span>
    <span class="k">const</span> <span class="n">MapVec</span> <span class="n">yy</span> <span class="o">=</span> <span class="n">Rcpp</span><span class="o">::</span><span class="n">as</span><span class="o"><</span><span class="n">MapVec</span><span class="o">></span><span class="p">(</span><span class="n">y</span><span class="p">);</span>
    <span class="c1">// Negative log likelihood
</span>    <span class="n">LogisticReg</span> <span class="n">nll</span><span class="p">(</span><span class="n">xx</span><span class="p">,</span> <span class="n">yy</span><span class="p">);</span>
    <span class="c1">// Initial guess
</span>    <span class="n">Eigen</span><span class="o">::</span><span class="n">VectorXd</span> <span class="n">beta</span><span class="p">(</span><span class="n">xx</span><span class="p">.</span><span class="n">cols</span><span class="p">());</span>
    <span class="n">beta</span><span class="p">.</span><span class="n">setZero</span><span class="p">();</span>

    <span class="kt">double</span> <span class="n">fopt</span><span class="p">;</span>
    <span class="kt">int</span> <span class="n">status</span> <span class="o">=</span> <span class="n">optim_lbfgs</span><span class="p">(</span><span class="n">nll</span><span class="p">,</span> <span class="n">beta</span><span class="p">,</span> <span class="n">fopt</span><span class="p">);</span>
    <span class="k">if</span><span class="p">(</span><span class="n">status</span> <span class="o"><</span> <span class="mi">0</span><span class="p">)</span>
        <span class="n">Rcpp</span><span class="o">::</span><span class="n">stop</span><span class="p">(</span><span class="s">"fail to converge"</span><span class="p">);</span>

    <span class="k">return</span> <span class="n">Rcpp</span><span class="o">::</span><span class="n">wrap</span><span class="p">(</span><span class="n">beta</span><span class="p">);</span>
<span class="p">}</span>

Now let’s do a quick benchmark:

<span class="n">set.seed</span><span class="p">(</span><span class="m">123</span><span class="p">)</span>
<span class="n">n</span> <span class="o">=</span> <span class="m">5000</span>
<span class="n">p</span> <span class="o">=</span> <span class="m">100</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">matrix</span><span class="p">(</span><span class="n">rnorm</span><span class="p">(</span><span class="n">n</span> <span class="o">*</span> <span class="n">p</span><span class="p">),</span> <span class="n">n</span><span class="p">)</span>
<span class="n">beta</span> <span class="o">=</span> <span class="n">runif</span><span class="p">(</span><span class="n">p</span><span class="p">)</span>
<span class="n">xb</span> <span class="o">=</span> <span class="n">c</span><span class="p">(</span><span class="n">x</span> <span class="o">%*%</span> <span class="n">beta</span><span class="p">)</span>
<span class="n">p</span> <span class="o">=</span> <span class="n">exp</span><span class="p">(</span><span class="n">xb</span><span class="p">)</span> <span class="o">/</span> <span class="p">(</span><span class="m">1</span> <span class="o">+</span> <span class="n">exp</span><span class="p">(</span><span class="n">xb</span><span class="p">))</span>
<span class="n">y</span> <span class="o">=</span> <span class="n">rbinom</span><span class="p">(</span><span class="n">n</span><span class="p">,</span> <span class="m">1</span><span class="p">,</span> <span class="n">p</span><span class="p">)</span>

<span class="n">system.time</span><span class="p">(</span><span class="n">res1</span> <span class="o"><-</span> <span class="n">glm.fit</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">family</span> <span class="o">=</span> <span class="n">binomial</span><span class="p">())</span><span class="o">$</span><span class="n">coefficients</span><span class="p">)</span>
<span class="c1">##  user  system elapsed
## 0.339   0.004   0.342
</span><span class="n">system.time</span><span class="p">(</span><span class="n">res2</span> <span class="o"><-</span> <span class="n">logistic_reg</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">))</span>
<span class="c1">##  user  system elapsed
##  0.01    0.00    0.01
</span><span class="n">max</span><span class="p">(</span><span class="n">abs</span><span class="p">(</span><span class="n">res1</span> <span class="o">-</span> <span class="n">res2</span><span class="p">))</span>
<span class="c1">## [1] 1.977189e-07
</span>

This is not a fair comparison however, since glm.fit() will calculate some
other components besides $beta$, and the precision of two methods are also
different.

RcppNumerical provides a function fastLR() that is a more stable
version of the code above (avoiding exp() overflow) and returns similar
components as glm.fit(). The performance is similar:

<span class="n">system.time</span><span class="p">(</span><span class="n">res3</span> <span class="o"><-</span> <span class="n">fastLR</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">)</span><span class="o">$</span><span class="n">coefficients</span><span class="p">)</span>
<span class="c1">##  user  system elapsed
##  0.01    0.00    0.01
</span><span class="n">max</span><span class="p">(</span><span class="n">abs</span><span class="p">(</span><span class="n">res1</span> <span class="o">-</span> <span class="n">res3</span><span class="p">))</span>
<span class="c1">## [1] 1.977189e-07
</span>

Its source code can be found
here.

Final Words

If you think this package may be helpful, feel free to leave comments or
request features in the Github
page. Contribution and pull requests would be great.

To leave a comment for the author, please follow the link and comment on their blog: Yixuan's Blog - R.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)