Rcpp only returning the first part of return line

Question

I'm attempting to code the the log-likelihood for e.g. linear mixed models seen in page 965 here.

I would implement this in R, quite trivially, as

R.imp <- function(Y, X, Z, V, b, beta, mi){
-mi/2 * log(2*pi) - 1/2 * log(det(V)) - 1/2 * crossprod(Y - X %*% beta - Z %*% b, solve(V) %*% (Y - X %*% beta - Z %*% b))
}

And my effort thus-far in implementing in Rcpp is

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace arma;
// [[Rcpp::export]]
double getLongit(vec& b, const colvec& Y, const mat& X, const mat& Z, const vec& beta, const mat& V, int mi){
    colvec resid = Y - X * beta - Z * b;
    // Rcpp::Rcout << "resid: " << resid << std::endl;
    // Rcpp::Rcout << "RTVR: " << resid.t() * V.i() * resid << std::endl;
    return as_scalar(-mi/2 * log(2 * M_PI) - 1/2 * log(det(V)) -1/2 * resid.t() * V.i() * resid);
}

where the commented lines are for debugging purposes (so clearly hasn't been too fruitful on my end).

A toy example to illustrate the problem:

set.seed(100)
Y <- rnorm(8); mi <- length(Y)
Z <- cbind(1, 0:7); b <- rnorm(2)
X <- cbind(1, rnorm(8), rbinom(8, 1, 0.5)); beta <- c(10, 5, -2)
V <- diag(8)

Which returns

> LongitR(b, Y, X, Z, beta, V, mi)
          [,1]
[1,] -232.7768
> getLongit(b, Y, X, Z, beta, V, mi)
[1] -7.351508

Where the result from my C++ implementation is the value of -8/2*log(2pi), the first term in the return statement and I'm unsure why nothing gets parsed beyond this?

I am assuming I've missed something quite obvious!

I recommend debugging it piece by piece. You have all the components, so check which one is falling short. — Dirk Eddelbuettel, Jun 24 '21 at 12:30
Your second term vanishes. I get zero for `as_scalar(-1/2 * resid.t() * V.i() * resid)`. — Dirk Eddelbuettel, Jun 24 '21 at 12:39
@DirkEddelbuettel thanks very much for checking - my Rcout lines didn't have the -1/2 term (looking ahead to your answer) and so these returned as they should without it, however lesson well and truly learnt! — jmurray, Jun 24 '21 at 16:34

Dirk Eddelbuettel · Accepted Answer · 2021-06-24T13:23:30.257

It is our least favourite error. You let integer math creep in via the 1/2 term!!

If you switch that to 0.5, all is well. A slightly edited version in a single file is below.

Output

> Rcpp::sourceCpp("~/git/stackoverflow/68115768/answer.cpp")

> getLL_R <- function(Y, X, Z, V, b, beta, mi) {
+     resid <- Y - X %*% beta - Z %*% b
+     -mi/2 * log(2*pi) - 1/2 * log(det(V)) - 1/2 * crossprod .... [TRUNCATED] 

> set.seed(100)

> Y <- rnorm(8)

> mi <- length(Y)

> Z <- cbind(1, 0:7)

> b <- rnorm(2)

> X <- cbind(1, rnorm(8), rbinom(8, 1, 0.5))

> beta <- c(10, 5, -2)

> V <- diag(8)

> getLL_R(Y, X, Z, V, b, beta, mi)
         [,1]
[1,] -232.777

> getLL_Cpp(b, Y, X, Z, beta, V, mi)
[1] -232.777
>

Code

#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]
using namespace arma;

// [[Rcpp::export]]
double getLL_Cpp(vec& b, const colvec& Y, const mat& X,
                 const mat& Z, const vec& beta, const mat& V, int mi){
    colvec resid = Y - X * beta - Z * b;
    return as_scalar(-mi/2.0 * log(2.0 * M_PI) -
                     0.5 * log(det(V)) - 0.5 * resid.t() * V.i() * resid);
}

/*** R
getLL_R <- function(Y, X, Z, V, b, beta, mi) {
    resid <- Y - X %*% beta - Z %*% b
    -mi/2 * log(2*pi) - 1/2 * log(det(V)) -
          1/2 * crossprod(resid, solve(V) %*% resid)
}

set.seed(100)
Y <- rnorm(8)
mi <- length(Y)
Z <- cbind(1, 0:7)
b <- rnorm(2)
X <- cbind(1, rnorm(8), rbinom(8, 1, 0.5))
beta <- c(10, 5, -2)
V <- diag(8)

getLL_R(Y, X, Z, V, b, beta, mi)
getLL_Cpp(b, Y, X, Z, beta, V, mi)
*/

D'oh - thanks so much. I've only really programmed in R up to this point, so the integer math creep is a new concept to me entirely. Lesson well and truly learnt, though - thanks a lot! — jmurray, Jun 24 '21 at 16:34
Trust me, _we all_ have wasted time on this. A rite of passage. Compilers got better about warning about this so it pays to always enable `-Wall -pedantic` and alike. — Dirk Eddelbuettel, Jun 24 '21 at 16:41
If I may, a minor bit of SO tradition: Thanks for accepting, please feel free to also upvote (click on up-arrow). Thanks! — Dirk Eddelbuettel, Jun 24 '21 at 16:42
I will certainly endeavour to once I reach the 15 reputation! — jmurray, Jun 24 '21 at 17:42
Sorry, I keep forgetting about that constraint. I just helped you along a little more as the question was actually well posted and framed and reproducible -- thumbs up! — Dirk Eddelbuettel, Jun 24 '21 at 17:47

Rcpp only returning the first part of return line

1 Answers1

Output

Code