Rcpp only returning the first part of return line

user_16306486 0
Viewed 48

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!

rcpparmadillo rcpp r
1 Answers

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)
*/
user_1433050
Related