У меня есть код C ++, выполняемый в Rcpp, где я определяю несколько функций, которые затем вызываются в открытой функции с помощью тега // [[Rcpp::export]].Код компилируется нормально, но выполнение открытой функции возвращает в фатальном сбое моего сеанса R, приводящем к немедленному завершению.
Что меня удивляет, так это то, что код выполнялся нормально вчера, когда я выполнил его до строки VectorXd z = y_luet - kroneckerProduct(X_luet.transpose(), MatrixXd::Identity(p, p)) * r;и возвращая вектор z.Теперь ни тот, ни полный код, как показано ниже, не работают.
Я также выполнил домашнее задание, протестировав все функции по отдельности, проверив их правильность, выставив им R перед использованием той же техники и проверив их по их более медленным аналогам R, получив численно идентичные результаты (при большемскорость).
Мне интересно, я просто не использую подход «определить несколько функций, а затем использовать их в большей функции», как только задачи становятся немного больше?
Сами данные являются умеренными по стандартам Эйгена, dat представляет собой матрицу с 200 строками и 2 столбцами, все остальное является низкоразмерным, с максимумом (строка, столбец) не более 12, т.е. вторая по величине матрица равна 12на 1.
Я использую Rtools и Rcpp для всех самых последних сборок.
В коде реализована простая итеративная обобщенная оценка наименьших квадратов, как это обычно делается в статистике / эконометрике.Я был бы очень признателен за вашу помощь.
РЕДАКТИРОВАТЬ: вот некоторые примеры данных в формате R, которые должны получить минимальный рабочий пример:
params <- .963
G <- matrix(c(1,0),nrow = 2)
G_perp <- matrix(c(0,1),nrow = 2)
mat_Lambda_lu <- matrix(0.95,nrow=1)
dat <- matrix(c(0,0,0,-0.79642284,-1.36694331,-1.18267593,-1.48827199,0.12549353,3.03343410,7.36256542,0,0,0,0.11282054,0.24798861,0.32448004,-0.27283699,-1.2462477,-0.0104694,3.21067339), nrow = 10, ncol = 2)
k<-3
no_ur <- 1
maxiter <- 100 #or something small for conserving memory
mini <- TRUE
Выше должно быть выполнено вR среда для его работы.Пожалуйста, дайте мне знать, если он не работает или есть проблемы.
Вот код:
// [[Rcpp::depends(RcppEigen)]]
#include <Rcpp.h>
#include <RcppEigen.h>
#include <cmath>
#include <cstdlib>
#include <Eigen/Dense>
#include <unsupported/Eigen/src/MatrixFunctions/MatrixPower.h>
#include <unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h>
using namespace Rcpp;
using namespace Eigen;
//using Eigen::Map; // 'maps' rather than copies
//using Eigen::MatrixXd; // variable size matrix, double precision
//using Eigen::VectorXd; // variable size vector, double precision
MatrixXd makeXluet(MatrixXd dat, int k, int p, int T) {
MatrixXd X_luet(p*k, T - k);
for (int i = k; i > 0; i--)
{
X_luet.block((k-i)*p,0,p,T-k) = dat.block(i - 1, 0, T - k, p).transpose();
}
return X_luet;
}
MatrixXd makeRLuTilde(MatrixXd Rlu, MatrixXd LambdaLu, int k, int p, int q) {
MatrixXd RLuTilde(p*k, q);
MatrixPower<MatrixXd> Apow(LambdaLu);
for (int i = k; i > 0; i--) {
RLuTilde.block((k - i)*p, 0, p, q) = Rlu * Apow(i-1);
}
return RLuTilde;
}
VectorXd GLSEstimateFast(MatrixXd Xluet, MatrixXd Sigma_u, MatrixXd R, VectorXd z, int T, int k) {
return (R.transpose() * kroneckerProduct(Xluet * Xluet.transpose(), Sigma_u.inverse()) * R).inverse() * R.transpose() * kroneckerProduct(Xluet, Sigma_u.inverse()) * z;
}
MatrixXd ResMaker(MatrixXd Xluet, MatrixXd Yluet, VectorXd beta, const int k, const int p) {
Map<MatrixXd> A(beta.data(), p, k*p);
return Yluet - A * Xluet;
}
double GLSCriterion(MatrixXd res, MatrixXd Sigma_u, const int k, const int p, const int T) {
MatrixXd Lp = Sigma_u.inverse().llt().matrixL().transpose();
MatrixXd v = Lp * res;
Map<VectorXd> v2(v.data(), v.size());
return (1 / static_cast<double>(T)) * v2.transpose() * v2;
}
MatrixXd CovEstFast(MatrixXd res, const int T) {
return (1 / static_cast<double>(T)) * res * res.transpose();
}
double likeli_h(MatrixXd CovEstHat, const int T) {
return (-0.5)*static_cast<double>(T) * log(CovEstHat.determinant());
}
// [[Rcpp::export]]
double restricted_iterated_ml_cpp(Map<VectorXd> params, Map<MatrixXd> G, Map<MatrixXd> G_perp, Map<MatrixXd> mat_Lambda_lu, Map<MatrixXd> dat, const int k, const int no_ur, const int maxiter, bool mini) {
const int p = dat.cols();
const int T = dat.rows();
int p2 = static_cast<int>(pow(p, 2));
int iter = 0;
MatrixXd X_luet = makeXluet(dat, k, p, T);
MatrixXd Y_luet = dat.bottomRows(T-k).transpose();
Map<MatrixXd> D(params.data(), p - no_ur, no_ur);
MatrixXd R_lu = G + G_perp * D;
MatrixXd R_lu_tilde = makeRLuTilde(R_lu, mat_Lambda_lu, k, p, no_ur);
MatrixXd C = kroneckerProduct(R_lu_tilde.transpose(), MatrixXd::Identity(T - k, T - k));
MatrixXd C1 = C.topLeftCorner(no_ur*p, no_ur*p);
MatrixXd C2 = C.block(0, no_ur*p, no_ur*p, C.cols() - (no_ur*p));
MatrixPower<MatrixXd> Llupow(mat_Lambda_lu);
MatrixXd mat_cee = R_lu * Llupow(k);
Map<VectorXd> cee(mat_cee.data(), mat_cee.size());
MatrixXd R(no_ur*p + k * p2 - (no_ur * p), k*p2 - (no_ur * p));
R << static_cast<double>(-1) * C1.inverse()*C2,
MatrixXd::Identity(k*p2-no_ur*p, k*p2 - (no_ur * p));
VectorXd r(k * p2);
r << C1.inverse() * cee,
MatrixXd::Zero(k * p2 - (no_ur * p), 1);
Map<VectorXd> y_luet(Y_luet.data(), Y_luet.size());
VectorXd z = y_luet - kroneckerProduct(X_luet.transpose(), MatrixXd::Identity(p, p)) * r;
MatrixXd Sigma_u = MatrixXd::Identity(p, p);
VectorXd gamma = GLSEstimateFast(X_luet, Sigma_u, R, z, T, k);
VectorXd beta = R * gamma + r;
MatrixXd res = ResMaker(X_luet, Y_luet, beta, k, p);
double crit_old = GLSCriterion(res, Sigma_u, k, p, T);
double crit_new = crit_old;
do
{
crit_old = crit_new;
Sigma_u = CovEstFast(res, T);
gamma = GLSEstimateFast(X_luet, Sigma_u, R, z, T, k);
beta = R * gamma + r;
res = ResMaker(X_luet, Y_luet, beta, k, p);
crit_new = GLSCriterion(res, Sigma_u, k, p, T);
iter++;
} while ((iter<maxiter) && (crit_old-crit_new>0.001));
double ll = likeli_h(Sigma_u, T);
if (mini) {
ll = static_cast<double>(-1)*ll;
}
return ll;
}