---
title: "Scores, Hessians, and varying parameters"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Scores, Hessians, and varying parameters}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(rvinecopulib)
set.seed(601)
```

rvinecopulib exposes observation-wise likelihood scores and average Hessians
for parametric bivariate and vine copula models. These are the basic inputs for
convergence checks, standard errors, sandwich covariance estimates, and
gradient-based extensions.

For a parameter vector `theta` and observation-wise log-likelihood
contributions `ell_i(theta)`, rvinecopulib reports
\[
s_i(\theta)=\frac{\partial\ell_i(\theta)}{\partial\theta},
\qquad
\bar H(\theta)=\frac{1}{n}\sum_{i=1}^n
\frac{\partial^2\ell_i(\theta)}{\partial\theta\partial\theta^\top}.
\]
Thus `scores()` contains the row vectors `s_i`, while `hessian()` returns
`H-bar`.

## Bivariate scores and Hessians

Fit a parametric model and evaluate derivatives at its fitted parameters:

```{r bicop-derivatives}
u2 <- rbicop(250, "gaussian", 0, 0.6)
fit2 <- bicop(u2, family_set = "gaussian")

score2 <- scores(u2, fit2)
hessian2 <- hessian(u2, fit2)
dim(score2)
hessian2
colSums(score2)
```

The score matrix has one row per observation and one column per model
parameter. At an interior maximum-likelihood estimate, its column sums should
be close to zero. The Hessian is averaged over observations; multiply by the
sample size when a summed log-likelihood Hessian is required.

These derivatives are available for continuous parametric models. Boundary
estimates and weakly identified parameters need the same care as in any
likelihood analysis.

## Vine parameter ordering

For a vine, score columns follow `(tree, edge, parameter)` order: tree 1 from
left to right, then tree 2, and so on; multi-parameter pair copulas contribute
their parameters in the order shown by `get_parameters()`.

```{r fit-vine}
n <- 220
z <- rnorm(n)
x <- cbind(
  z + rnorm(n),
  0.7 * z + rnorm(n),
  -0.5 * z + rnorm(n)
)
u <- pseudo_obs(x)
fit <- vinecop(
  u,
  structure = dvine_structure(1:3),
  family_set = "gaussian"
)
get_all_parameters(fit)
```

The number and ordering of derivative columns can be checked against the edge
summary.

```{r vine-scores}
score <- scores(u, fit)
hess <- hessian(u, fit)
dim(score)
dim(hess)
summary(fit)
```

## Stepwise and full likelihood derivatives

Vine models are usually fitted sequentially. With `step_wise = TRUE` (the
default), each pair copula treats the pseudo-observations passed down from
earlier trees as fixed. This is the objective optimized by the sequential
estimator.

With `step_wise = FALSE`, derivatives propagate through the complete
h-function cascade and refer to the full joint log-likelihood.

For a vine with edge set `E`, an observation contributes
\[
\ell_i(\theta)=\sum_{e\in E}
\log c_e\{u_{e,1,i}(\theta_{<e}),u_{e,2,i}(\theta_{<e});\theta_e\}.
\]
Stepwise derivatives hold the conditional pseudo-observations `u_e` fixed.
Full derivatives also differentiate their dependence on parameters in earlier
trees, denoted by `theta_<e` above.

```{r stepwise-full}
stepwise_gradient <- colSums(scores(u, fit, step_wise = TRUE))
full_gradient <- colSums(scores(u, fit, step_wise = FALSE))
rbind(stepwise = stepwise_gradient, full = full_gradient)
```

The stepwise gradient should be close to zero at a sequential fit. The full
gradient generally is not: a sequence of conditional maximizations is not the
same as a joint maximum. This distinction should be stated explicitly whenever
derivatives are used for inference.

## Reuse density intermediates

`dvinecop(..., keep_all = TRUE)` returns the density and the per-edge values
computed by the h-function cascade.

```{r keep-all}
full <- dvinecop(u[1:5, ], fit, keep_all = TRUE)
names(full)
full$pdf
```

The triangular `pdf_edges`, `hfunc1`, and `hfunc2` entries follow the same tree
and edge ordering as the fitted pair copulas. The `_sub` h-functions contain
left-limit evaluations for discrete models and are empty for fully continuous
models.

## Observation-specific bivariate parameters

Parametric copulas can be evaluated with a different parameter set in every
row. This is useful when another model maps covariates to valid copula
parameters.

```{r varying-bicop}
points <- matrix(runif(200), ncol = 2)
rho <- seq(-0.8, 0.8, length.out = nrow(points))

density <- dbicop(points, "gaussian", 0, parameters = rho)
score_varying <- scores(points, bicop_dist("gaussian"), parameters = rho)
hessian_varying <- hessian(
  points,
  bicop_dist("gaussian"),
  parameters = rho
)
head(cbind(rho, density, score_varying))
```

For a one-parameter family, `parameters` may be a vector. Otherwise it must be
a matrix with one row per observation and one column per family parameter.
Parameters are not recycled.

## Observation-specific vine parameters

The vine interface takes an `n` by `p` parameter matrix, where `p` is the total
number of pair-copula parameters and columns use the score ordering.

```{r varying-vine}
base_parameters <- unlist(get_all_parameters(fit), use.names = FALSE)
parameter_matrix <- matrix(
  rep(base_parameters, each = nrow(u)),
  nrow = nrow(u)
)
parameter_matrix[, 1] <- seq(-0.6, 0.6, length.out = nrow(u))

density_varying <- dvinecop(u, fit, parameters = parameter_matrix)
score_varying <- scores(u, fit, parameters = parameter_matrix)
head(density_varying)
dim(score_varying)
```

Observation-specific vine parameters are supported for continuous,
all-parametric models. Each row must represent a valid complete parameter
vector. The model stored in `fit` is not modified.

## From derivatives to uncertainty estimates

The score and Hessian outputs deliberately provide ingredients rather than a
single universal covariance estimator. The correct assembly depends on whether
the estimator is stepwise or joint, whether weights or clusters are present,
and whether robust inference is desired.

For an ordinary interior bivariate MLE, the inverse negative summed Hessian is
the familiar model-based covariance estimate. More elaborate vine inference
should account for the sequential estimation scheme instead of treating the
stepwise estimate as a joint maximum.

In the ordinary independent-observation case, define
\[
A=-\bar H, \qquad B=\frac{1}{n}\sum_{i=1}^n s_i s_i^\top.
\]
The model-based covariance estimate is `(-n H-bar)^(-1)`, while the usual
sandwich estimate is `A^(-1) B A^(-1) / n`. Weights, clusters, and sequential
estimation change how these pieces should be assembled.

## Related documentation

- [Bivariate copula models](bivariate-copulas.html) and [vine copula models and
  structures](vine-copula-models.html) provide the model definitions used in
  the examples.
- [Conditional simulation and Rosenblatt
  transforms](conditional-simulation.html) explains the h-function cascade
  that full derivatives propagate through.
- The [bivariate distribution reference](../reference/bicop_methods.html) and
  [vine-copula distribution reference](../reference/vinecop_methods.html)
  specify derivative support and parameter shapes.