Rcpp 仅 returning return 行的第一部分

Question

我正在尝试编写对数似然的代码，例如page 965 here.

中的线性混合模型

我会在 R 中实现这个，非常简单，因为

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))
}

到目前为止，我在 Rcpp 中实现的努力是

#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);
}

其中注释行用于调试目的（所以显然对我而言并没有太大成果）。

一个玩具示例来说明问题：

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)

哪个returns

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

我的 C++ 实现的结果是 -8/2*log(2pi) 的值，即 return 语句中的第一项，我不确定为什么除了这个之外什么都没有被解析？

我假设我错过了一些非常明显的东西！

Answer 1

这是我们最不喜欢的错误。你让整数数学通过 1/2 术语悄悄进入！！

如果将其切换为 0.5，一切都很好。下面是单个文件中略微编辑的版本。

输出

> Rcpp::sourceCpp("~/git/Whosebug/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
>

代码

#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)
*/

Rcpp 仅 returning return 行的第一部分

Rcpp only returning the first part of return line

r

rcpp

rcpparmadillo

输出

代码