我已尝试在R 中实现一个简单的隐马尔可夫模型（Hidden Markov Modell），不同状态使用不同的分布。我尝试使用TMB包来实现带有随机效应的C++实现。我已经发布了我的TMB代码和一个在R 中的示例。在MakeADFun() 之后我的代码会产生以下错误信息：

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"

R会话被中止。R遇到了严重错误。会话被终止。当我在模型中加入样条函数时，这种情况就会发生，所以我认为TMB代码中的GMRF() 可能有问题。如果完全不包含随机效应，MakeADFun() 就能工作。

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

解决方案

问题似乎在于你使用 S_emis 来创建 beta_sigma_re 的精度矩阵，尽管它本应只用于 beta_mu_re 的精度矩阵。

如果我理解正确，你的模型中 beta_mu_re 是对一个带惩罚项的样条的基底系数，用来建模 Y 与 time 之间的关系，它们服从分布 beta_mu_re ~ MVN(0, (lambda1 S_emis)-1)，其中lambda1是平滑参数，S_emis 是由mgcv给出的一个精度矩阵。基底维度（这里由mgcv自动选择）为9，因此 beta_mu_re 的长度为9，S_emis 是一个9×9的矩阵。

但你把 beta_sigma_re 定义为长度为1 的向量，因为方差参数没有随机效应。当你写下

        penalty += Type(0.5) * exp(log_lambda_emis(1)) *
            density::GMRF(S_emis).Quadform(beta_sigma_re.col(k));

你是在假设一个形如 beta_sigma_re ~ MVN(0, (lambda2 S_emis)-1) 的模型，但这行不通，因为长度为1 的向量不可能有一个9×9的精度矩阵。

如何解决这个问题取决于你打算如何使用这段代码。

• 作为起点，你可以检查，当上面两行代码被注释掉时，MakeADFun() 可以运行。
• 为了使代码具有通用性和灵活性，你需要为不同的参数存储不同的精度矩阵（例如，在你的示例中用于 mu 和 sigma 的矩阵）。
• 当没有随机效应时该怎么办？我确信有多种做法。在hmmTMB中，只有在存在随机效应时才计算惩罚项（Github链接）。