Rcpp 只返回返回行的第一部分 输出代码

如何解决Rcpp 只返回返回行的第一部分 输出代码

我正在尝试编码对数似然，例如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)

哪个返回

> LongitR(b,Y,mi)
          [,1]
[1,] -232.7768
> getLongit(b,mi)
[1] -7.351508

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

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

解决方法

这是我们最不喜欢的错误。你让整数数学通过 1/2 项渗透进来！！

如果您将其切换为 0.5，则一切正常。下面是单个文件中稍作修改的版本。

输出

> 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,mi)
         [,1]
[1,] -232.777

> getLL_Cpp(b,Y,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,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,0.5))
beta <- c(10,-2)
V <- diag(8)

getLL_R(Y,mi)
getLL_Cpp(b,mi)
*/