Scores, gradients and Hessians

Bivariate copulas

• scores(u) — an matrix of per-observation scores, .
• gradient(u) — the summed score, i.e. the gradient of the log-likelihood.
• hessian(u) — the Hessian of the log-likelihood.
• hessian_full(u) — the per-observation Hessians, not summed.
• scores_cov(u) — , the outer-product ("BHHH") covariance of the scores.
• scores_full(u) — the scores together with the intermediate quantities.

At a maximum-likelihood fit the gradient is numerically zero, which is the cheapest available check that a fit converged.

  Bicop model(
    BicopFamily::student, 0, (Eigen::VectorXd(2) << 0.5, 4.0).finished());
  auto u = model.simulate(500, false, { 1 });

  // n x npars matrix of per-observation scores: d log c / d theta
  Eigen::MatrixXd s = model.scores(u);

  // the summed score, i.e. the gradient of the log-likelihood
  Eigen::VectorXd g = model.gradient(u);

  // npars x npars Hessian of the log-likelihood, and the per-observation
  // outer-product covariance of the scores
  Eigen::MatrixXd h = model.hessian(u);
  Eigen::MatrixXd cov = model.scores_cov(u);

  std::cout << "npars: " << s.cols() << ", gradient norm: " << g.norm()
            << std::endl;

  // At the MLE the gradient is (numerically) zero, which is the cheapest
  // available check that a fit converged.
  Bicop fitted(BicopFamily::student);
  fitted.fit(u);
  std::cout << "gradient norm at the MLE: " << fitted.gradient(u).norm()
            << std::endl;

  // A sandwich standard error: (-H)^-1 cov (-H)^-1
  Eigen::MatrixXd hinv = (-h).inverse();
  Eigen::MatrixXd sandwich = hinv * cov * hinv;
  std::cout << "std. errors: " << sandwich.diagonal().cwiseSqrt().transpose()
            << std::endl;

Vine copulas

The vine versions take the same names. The columns of scores(), and hence the entries of gradient() and the rows and columns of hessian(), are ordered (tree, edge, parameter): tree 0's edges from left to right, and within an edge the pair copula's own parameters in order.

The step_wise flag selects which derivative you get:

• step_wise = true (the default) gives the score of the step-wise MLE. Each pair copula's gradient treats its pseudo-observations as fixed, which is what the sequential fitting procedure actually optimizes.
• step_wise = false gives the gradient of the full log-likelihood, propagating through the h-function cascade.

The distinction is not cosmetic: at a fitted model the step-wise gradient vanishes while the full gradient generally does not, because the sequential procedure is not a joint maximizer.

  auto data = dependent_data(500, 4, 2);
  FitControlsVinecop controls(bicop_families::parametric);
  Vinecop model(data, RVineStructure(), {}, controls);

  // Columns are ordered (tree, edge, parameter): tree 0's edges first, and
  // within an edge the pair copula's own parameters in order.
  Eigen::MatrixXd s = model.scores(data);
  Eigen::VectorXd g = model.gradient(data);
  Eigen::MatrixXd h = model.hessian(data);

  std::cout << "total parameters: " << s.cols()
            << ", gradient norm: " << g.norm() << std::endl;

  // step_wise = true (the default) is the score of the step-wise MLE: each
  // pair copula's gradient treats its pseudo-observations as fixed. Pass false
  // for the gradient of the full log-likelihood, which propagates through the
  // h-function cascade.
  Eigen::VectorXd g_full = model.gradient(data, false);
  std::cout << "step-wise norm: " << g.norm()
            << ", full norm: " << g_full.norm() << std::endl;

Reusing intermediate quantities

The derivative cascade is the expensive part. scores_full() returns the per-edge densities, h-functions and their derivatives alongside the scores, and hessian_full() returns the per-edge Hessian blocks rather than assembling them, so a caller that needs several of these pays for the cascade once.

  auto data = dependent_data(300, 4, 3);
  FitControlsVinecop controls(bicop_families::parametric);
  Vinecop model(data, RVineStructure(), {}, controls);

  // scores_full() returns the intermediate per-edge quantities alongside the
  // scores, so a caller needing several of them pays the cascade once.
  auto full = model.scores_full(data);
  std::cout << "scores: " << full.scores.rows() << " x " << full.scores.cols()
            << std::endl;

  // hessian_full() keeps the per-edge blocks rather than assembling them
  auto blocks = model.hessian_full(data);

Per-observation parameters

pdf(), pdf_full(), scores(), scores_full(), gradient() and hessian() also accept an matrix holding one full-vine parameter vector per observation, with columns in the same (tree, edge, parameter) order. Broadcasting a single vector across the rows reproduces the fixed-parameter results exactly; varying the rows is what makes covariate-dependent parameters possible. Continuous, all-parametric models only.

  auto data = dependent_data(200, 3, 4);
  FitControlsVinecop controls(bicop_families::parametric);
  Vinecop model(data, RVineStructure(), {}, controls);

  // Flatten the stored parameters in the same (tree, edge, parameter) order
  // that scores() reports and the `parameters` argument expects.
  std::vector<double> flat;
  for (const auto& tree : model.get_all_parameters()) {
    for (const auto& edge : tree) {
      for (Eigen::Index p = 0; p < edge.size(); ++p) {
        flat.push_back(edge(p));
      }
    }
  }
  Eigen::VectorXd theta = Eigen::Map<Eigen::VectorXd>(
    flat.data(), static_cast<Eigen::Index>(flat.size()));

  // One parameter vector per observation. Broadcasting the stored parameters
  // reproduces the fixed-parameter results; varying the rows is what makes
  // covariate-dependent parameters possible. Continuous, all-parametric models
  // only.
  Eigen::MatrixXd parameters = theta.transpose().replicate(data.rows(), 1);

  auto pdf = model.pdf(data, parameters);
  auto s = model.scores(data, parameters);
  std::cout << "parameters per row: " << parameters.cols()
            << ", pdf matches the fixed-parameter path: "
            << pdf.isApprox(model.pdf(data)) << std::endl;