ARTICLE AD BOX
I have tried to implement a simple Hidden Markov Modell with different distributions for different states in R. I tried to use the TMB package to implement random effects with C++. I have posted my TMB code and an example in R. After MakeADFun() my code produces the following error messsage:
TMB has received an error from Eigen. The following condition was not met: lhs.cols()==rhs.cols() && "invalid mmatrix product" && "if you wanted coefficientwise or a dot product use the respective explicit functions"with R Session Aborted. R encountered a fatal error. The session was terminated. This happens as soon as I include splines in my modell, so i think there could be something wrong with GMRF() in my TMB code. If I include no effects at all MakeADFun() works.
#include <TMB.hpp> // =============================== // Log-density per state // =============================== template<class Type> Type log_density(Type x, Type mu, Type sigma, int dist_code) { if(dist_code == 0){ // Normal return dnorm(x, mu, sigma, true); } else if(dist_code == 1){ // Lognormal return dnorm(log(x), mu, sigma, true) - log(x); } else if(dist_code == 2){ // Gamma Type shape = pow(mu / sigma, 2); Type scale = pow(sigma, 2) / mu; return dgamma(x, shape, scale, true); } else if(dist_code == 3){ // Weibull return dweibull(x, sigma, mu, true); // shape = sigma, scale = mu } return Type(0); } // =============================== // HMM with Splines + Fixed Effects // =============================== template<class Type> Type objective_function<Type>::operator() () { // ================= // DATA // ================= DATA_VECTOR(y); DATA_INTEGER(K); DATA_IVECTOR(dist); int T = y.size(); // --- Emission Design --- DATA_SPARSE_MATRIX(X_mu_fe); DATA_SPARSE_MATRIX(X_mu_re); DATA_SPARSE_MATRIX(X_sigma_fe); DATA_SPARSE_MATRIX(X_sigma_re); // --- TPM Design --- DATA_SPARSE_MATRIX(X_tpm_fe); DATA_SPARSE_MATRIX(X_tpm_re); // --- Penalty Matrices --- DATA_SPARSE_MATRIX(S_emis); DATA_SPARSE_MATRIX(S_tpm); DATA_VECTOR(log_lambda_emis); DATA_VECTOR(log_lambda_tpm); // ================= // PARAMETERS // ================= // Emissions starting values PARAMETER_VECTOR(mu0); PARAMETER_VECTOR(sigma0); // Emission effects PARAMETER_MATRIX(beta_mu_fe); PARAMETER_MATRIX(beta_mu_re); PARAMETER_MATRIX(beta_sigma_fe); PARAMETER_MATRIX(beta_sigma_re); // Transition effects PARAMETER_MATRIX(beta_tpm_fe); PARAMETER_MATRIX(beta_tpm_re); // initialdistribution PARAMETER_VECTOR(delta_raw); // ================= // initialdistribution // ================= vector<Type> delta = exp(delta_raw); delta /= delta.sum(); // ================= // EMISSIONS // ================= matrix<Type> mu(T,K); matrix<Type> sigma(T,K); for(int t=0; t<T; t++){ for(int k=0; k<K; k++){ mu(t,k) = mu0(k) + (X_mu_fe.row(t) * beta_mu_fe.col(k)).sum() + (X_mu_re.row(t) * beta_mu_re.col(k)).sum(); sigma(t,k) = exp( log(sigma0(k)) + (X_sigma_fe.row(t) * beta_sigma_fe.col(k)).sum() + (X_sigma_re.row(t) * beta_sigma_re.col(k)).sum() ); } } // ================= // TRANSITIONS // ================= array<Type> Gamma(T-1,K,K); for(int t=0; t<T-1; t++){ for(int i=0; i<K; i++){ vector<Type> eta(K); for(int j=0; j<K; j++){ int idx = i*K + j; eta(j) = (X_tpm_fe.row(t) * beta_tpm_fe.col(idx)).sum() + (X_tpm_re.row(t) * beta_tpm_re.col(idx)).sum(); } Type max_eta = eta.maxCoeff(); vector<Type> exp_eta = (eta - max_eta).array().exp(); Type denom = exp_eta.sum(); for(int j=0; j<K; j++){ Gamma(t,i,j) = exp_eta(j) / denom; } } } // ================= // LOG-EMISSIONS // ================= matrix<Type> log_emis(T,K); for(int t=0; t<T; t++){ for(int k=0; k<K; k++){ log_emis(t,k) = log_density(y(t), mu(t,k), sigma(t,k), dist(k)); } } // ================= // FORWARD ALGORITHM // ================= matrix<Type> alpha(T,K); vector<Type> tmp(K); // t = 0 for(int k=0; k<K; k++){ tmp(k) = delta(k) * exp(log_emis(0,k)); } Type scale = tmp.sum(); if(scale < Type(1e-16)) scale = Type(1e-16); alpha.row(0) = (tmp / scale).transpose(); Type logLik = log(scale); // recursion for(int t=1; t<T; t++){ for(int j=0; j<K; j++){ Type s = 0; for(int i=0; i<K; i++){ s += alpha(t-1,i) * Gamma(t-1,i,j); } tmp(j) = s * exp(log_emis(t,j)); } scale = tmp.sum(); if(scale < Type(1e-16)) scale = Type(1e-16); alpha.row(t) = (tmp / scale).transpose(); logLik += log(scale); } // ================= // PENALTIES // ================= Type penalty = 0; // Emission Splines for(int k=0; k<K; k++){ penalty += Type(0.5) * exp(log_lambda_emis(0)) * density::GMRF(S_emis).Quadform(beta_mu_re.col(k)); penalty += Type(0.5) * exp(log_lambda_emis(1)) * density::GMRF(S_emis).Quadform(beta_sigma_re.col(k)); } // TPM Splines for(int j=0; j<K*K; j++){ penalty += Type(0.5) * exp(log_lambda_tpm(0)) * density::GMRF(S_tpm).Quadform(beta_tpm_re.col(j)); } // ================= // OBJECTIVE // ================= Type nll = -logLik + penalty; // ================= // REPORT // ================= REPORT(mu); REPORT(sigma); REPORT(Gamma); REPORT(alpha); ADREPORT(delta); return nll; } compile("hmm_full.cpp") dyn.load(TMB::dynlib("hmm_full")) library(hmmTMB) library(Matrix) Y <- rnorm(48) # Beispiel-Daten time <- 1:48 K <- 2 dist <- c(0,0) # beide Normal T <- length(Y) #no fixed effects for mu X_mu_fe <- Matrix(0, nrow=T, ncol=1, sparse=TRUE) #model mu as spline dependent on time ts_setup = gam(Y ~ s(time, bs = "ts") - 1, data = df, fit = FALSE) S = ts_setup$smooth[[1]]$S[[1]] X = ts_setup$X X_mu_re <- X # Designmatrix S_emis <- S # Penalty X_mu_re <- as_sparse(X_mu_re) S_emis <- as_sparse(S) #no fixed effects for sigma or tpm X_sigma_fe <- Matrix(0, nrow=T, ncol=1, sparse=TRUE) X_tpm_fe <- Matrix(0, nrow=T, ncol=1, sparse=TRUE) # --- no spline effects for sigma or tpm --- X_sigma_re <- Matrix(0, nrow=T, ncol=1, sparse=TRUE) X_tpm_re <- Matrix(0, nrow=T, ncol=1, sparse=TRUE) # --- dummy penalty matrix for tpm (no effects for tpm) --- S_tpm <- Matrix(1,1,1,sparse=TRUE) # --- starting values --- mu0 <- c(0,1) sigma0 <- c(0.5,0.8) beta_mu_fe <- matrix(0, ncol(X_mu_fe), K) beta_mu_re <- matrix(0, ncol(X_mu_re), K) beta_sigma_fe <- matrix(0, ncol(X_sigma_fe), K) beta_sigma_re <- matrix(0, ncol(X_sigma_re), K) beta_tpm_fe <- matrix(0, ncol(X_tpm_fe), K*K) beta_tpm_re <- matrix(0, ncol(X_tpm_re), K*K) delta_raw <- rep(0,K) log_lambda_emis <- c(0,0) log_lambda_tpm <- 0 data <- list( y=Y, K=K, dist=dist, X_mu_fe=X_mu_fe, X_mu_re=X_mu_re, X_sigma_fe=X_sigma_fe, X_sigma_re=X_sigma_re, X_tpm_fe=X_tpm_fe, X_tpm_re=X_tpm_re, S_emis=S_emis, S_tpm=S_tpm, log_lambda_emis=log_lambda_emis, log_lambda_tpm=log_lambda_tpm ) parameters <- list( mu0=mu0, sigma0=sigma0, beta_mu_fe=beta_mu_fe, beta_mu_re=beta_mu_re, beta_sigma_fe=beta_sigma_fe, beta_sigma_re=beta_sigma_re, beta_tpm_fe=beta_tpm_fe, beta_tpm_re=beta_tpm_re, delta_raw=delta_raw ) obj <- MakeADFun(data, parameters, DLL="hmm_full", silent=TRUE)