1 The pure inverse problem

Suppose we face the following inverse problem

\[\begin{equation} \mathbf{y} = \mathbf{X}\mathbf{p} \tag{1.1} \end{equation}\]

where \(\mathbf{y} = (y_1, \dots, y_T)'\) is a \(T\)-dimensional vector of observations (moments), \(\mathbf{X}\) is a known \((T \times K)\) matrix, and \(\mathbf{p} = (p_1, \dots, p_K)'\) is a \(K\)-dimensional vector of unknown probabilities. The goal is to recover \(\mathbf{p}\), but in many cases this task is not directly feasible — for example, when \(\mathbf{X}\) has no inverse. A classical example is the under-determined case in which the number of observations is smaller than the number of unknowns, \(T < K\). There are then infinitely many \(\mathbf{p}\) consistent with (1.1), so we need a criterion to select among them.

Throughout, inverse_ce() solves exactly this problem through its concentrated (dual) formulation.

library(infometrics)

2 Maximum entropy

Following Shannon (1948), (Jaynes, 1957a, 1957b), and Levine (1980), one natural criterion is to choose the \(\mathbf{p}\) that maximises the Shannon entropy

\[\begin{equation} H(\mathbf{p}) = -\sum_{k=1}^K p_k \log(p_k), \tag{2.1} \end{equation}\]

with the convention \(p_k \log(p_k) = 0\) whenever \(p_k = 0\). Entropy (2.1) is maximised when \(\mathbf{p}\) is uniform, i.e. \(p_k = 1/K\) for all \(k\), so maximising it selects the least committal distribution consistent with the data. We also impose the additivity (normalization) constraint

\[\begin{equation} \sum_{k=1}^K p_k = 1. \tag{2.2} \end{equation}\]

We therefore maximise (2.1) subject to (1.1) and (2.2). Writing the moment residual as \(\big(\sum_k p_k x_{tk} - y_t\big)\), the Lagrangian is

\[\begin{equation} \mathcal{L} = -\sum_{k=1}^K p_k \log(p_k) + \sum_{t=1}^T \lambda_t\!\left(\sum_{k=1}^K p_k x_{tk} - y_t\right) + \mu\!\left(1 - \sum_{k=1}^K p_k\right), \tag{2.3} \end{equation}\]

where the \(\lambda_t\) and \(\mu\) are Lagrange multipliers. The multiplier \(\lambda_t\) measures the information content of moment \(t\): if \(\lambda_t = 0\), the \(t\)-th moment carries no information and dropping it would not change the solution or the level of remaining uncertainty. The first-order conditions are

\[\begin{equation} \begin{aligned} \frac{\partial \mathcal{L}}{\partial p_k} &= -\log(p_k) - 1 + \sum_{t=1}^T \lambda_t x_{tk} - \mu = 0, & k &= 1,\dots,K,\\ \frac{\partial \mathcal{L}}{\partial \lambda_t} &= \sum_{k=1}^K p_k x_{tk} - y_t = 0, & t &= 1,\dots,T,\\ \frac{\partial \mathcal{L}}{\partial \mu} &= 1 - \sum_{k=1}^K p_k = 0. \end{aligned} \tag{2.4} \end{equation}\]

Solving the first condition and imposing (2.2) gives the exponential-family solution

\[\begin{equation} \hat{p}_k = \frac{\exp\!\big(\sum_t \hat\lambda_t x_{tk}\big)} {\sum_{k'} \exp\!\big(\sum_t \hat\lambda_t x_{tk'}\big)} = \frac{\exp\!\big(\sum_t \hat\lambda_t x_{tk}\big)}{\Omega(\hat{\boldsymbol\lambda})}, \tag{2.5} \end{equation}\]

where \(\Omega(\boldsymbol\lambda) = \sum_k \exp\!\big(\sum_t \lambda_t x_{tk}\big)\) is the normalization factor, usually called the partition function. Substituting (2.5) back into (2.3) — the normalization is already satisfied, so the \(\mu\) term drops — yields the concentrated (dual) model

\[\begin{equation} M(\boldsymbol\lambda) = -\sum_t \lambda_t y_t + \log \Omega(\boldsymbol\lambda), \tag{2.6} \end{equation}\]

which is minimised over the \(T\)-dimensional vector \(\boldsymbol\lambda\). This is exactly the objective inverse_ce() solves, so the reported lambda_hat and objective correspond one-to-one with \(\hat{\boldsymbol\lambda}\) and \(M(\hat{\boldsymbol\lambda})\) here. Working in \(\boldsymbol\lambda\)-space (dimension \(T\)) rather than \(\mathbf{p}\)-space (dimension \(K \gg T\)) is the key computational advantage of the dual (Golan, 2008).

dice <- data.frame(y = 4.5,
                   s1 = 1, s2 = 2, s3 = 3, s4 = 4, s5 = 5, s6 = 6)
fit_me <- inverse_ce(y ~ s1 + s2 + s3 + s4 + s5 + s6 - 1, data = dice)
round(coef(fit_me), 4)          # p_hat over the six faces
#>     s1     s2     s3     s4     s5     s6 
#> 0.0544 0.0788 0.1142 0.1654 0.2398 0.3475

sum(1:6 * coef(fit_me))          # reproduces the mean 4.5
#> [1] 4.5
c(H_phat = shannon_entropy(coef(fit_me)), H_max = log(6))
#>   H_phat    H_max 
#> 1.613581 1.791759

3 Cross-entropy: incorporating a prior

If we hold prior beliefs about the data-generating process, summarised by a distribution \(\mathbf{p}^0 = (p_1^0, \dots, p_K^0)'\), we can incorporate them through the cross-entropy (Kullback–Leibler) divergence (Csiszar, 1991; Kullback, 1974; Levine, 1980; Shore & Johnson, 1980)

\[\begin{equation} D(\mathbf{p} \,\Vert\, \mathbf{p}^0) = \sum_{k=1}^K p_k \log(p_k / p_k^0), \tag{3.1} \end{equation}\]

which measures the divergence of \(\mathbf{p}\) from the prior \(\mathbf{p}^0\). We now minimise (3.1) subject to the same constraints (1.1) and (2.2). The analogous derivation gives

\[\begin{equation} \hat{p}_k = \frac{p_k^0 \exp\!\big(\sum_t \hat\lambda_t x_{tk}\big)} {\sum_{k'} p_{k'}^0 \exp\!\big(\sum_t \hat\lambda_t x_{tk'}\big)} = \frac{p_k^0 \exp\!\big(\sum_t \hat\lambda_t x_{tk}\big)}{\Omega(\hat{\boldsymbol\lambda})}, \tag{3.2} \end{equation}\]

with \(\Omega(\boldsymbol\lambda) = \sum_k p_k^0 \exp\!\big(\sum_t \lambda_t x_{tk}\big)\). Maximum entropy is the special case of (3.2) with a uniform prior \(p_k^0 = 1/K\): the constant \(1/K\) cancels between numerator and denominator and (3.2) collapses to (2.5). The concentrated dual has the same form as (2.6), only the partition function changes, so inverse_ce() needs a single implementation for both: uniform p0 gives ME, a non-uniform p0 gives CE (hence the name).

Since the divergence (3.1) is non-negative, a strictly positive value means the data carried information beyond the prior. Following (Zellner, 1988, 1991), (3.2) can be read as an efficient information-processing rule that combines prior and data into a posterior.

fit_unif <- inverse_ce(y ~ s1 + s2 + s3 + s4 + s5 + s6 - 1,
                       data = dice, p0 = rep(1 / 6, 6))
max(abs(coef(fit_unif) - coef(fit_me)))     # ~ 0: ME == CE(uniform prior)
#> [1] 0

p0_load <- c(.05, .05, .10, .15, .25, .40)  # prior beliefs favouring high faces
fit_ce  <- inverse_ce(y ~ s1 + s2 + s3 + s4 + s5 + s6 - 1,
                      data = dice, p0 = p0_load)
round(rbind(ME = coef(fit_me), CE = coef(fit_ce)), 4)
#>        s1     s2     s3     s4     s5     s6
#> ME 0.0544 0.0788 0.1142 0.1654 0.2398 0.3475
#> CE 0.0689 0.0630 0.1153 0.1582 0.2413 0.3533
sum(1:6 * coef(fit_ce))                      # still satisfies the mean moment
#> [1] 4.5

4 The information matrix and standard errors

The gradient of the dual (2.6) returns the moment condition (1.1), and its Hessian is the Fisher information matrix for \(\boldsymbol\lambda\),

\[\begin{equation} \big[\mathbf{I}(\boldsymbol\lambda)\big]_{st} = \frac{\partial^2 M}{\partial \lambda_s\, \partial \lambda_t} = \sum_k x_{sk} x_{tk}\, \hat p_k - \Big(\sum_k x_{sk}\hat p_k\Big)\Big(\sum_k x_{tk}\hat p_k\Big) = \mathrm{Cov}_{\hat p}(\mathbf{x}_s, \mathbf{x}_t), \tag{4.1} \end{equation}\]

i.e. the covariance of the moment functions under the fitted distribution \(\hat{\mathbf p}\) (Golan Eq. 4.7). The diagonal elements are the variances \(\mathrm{Var}_{\hat p}(\mathbf{x}_t)\); the off-diagonal elements (\(s \neq t\)) are the covariances \(\mathrm{Cov}_{\hat p}(\mathbf{x}_s, \mathbf{x}_t)\). As shown in Chapter 17 of Cover & Thomas (2006) and Golan (2008) (p. 59), the same matrix can be recovered from the local curvature of the cross-entropy. The variance–covariance matrix of \(\hat{\boldsymbol\lambda}\) is its inverse,

\[\begin{equation} \mathrm{Var}(\hat{\boldsymbol\lambda}) = \mathbf{I}^{-1}(\boldsymbol\lambda), \tag{4.2} \end{equation}\]

and a delta-method step propagates it to \(\hat{\mathbf p}\). These are curvature / identification quantities, not sampling standard errors: a pure inverse problem is deterministic, so they describe how sharply the moments pin down \(\boldsymbol\lambda\), not sampling variability. For sampling inference, use the stochastic-moment sibling inverse_noise().

faces  <- 1:6
p_ref  <- c(.10, .12, .15, .18, .20, .25)
Xm     <- rbind(faces, faces^2)                 # 2 moments x 6 states
ym     <- as.numeric(Xm %*% p_ref)              # feasible (mean, 2nd moment)
d2 <- data.frame(y  = ym,
                 s1 = Xm[, 1], s2 = Xm[, 2], s3 = Xm[, 3],
                 s4 = Xm[, 4], s5 = Xm[, 5], s6 = Xm[, 6])
fit2 <- inverse_ce(y ~ s1 + s2 + s3 + s4 + s5 + s6 - 1, data = d2)
summary(fit2)
#> 
#> Call:
#> inverse_ce(formula = y ~ s1 + s2 + s3 + s4 + s5 + s6 - 1, data = d2)
#> 
#> Moment-fitting residuals (y - Xp):
#> [1]  -3.888e-07  -3.081e-06
#> 
#> Coefficients (estimated probabilities p; delta-method SE from I^-1):
#>     Estimate  Std. Error
#> s1    0.1002      0.2583
#> s2    0.1218      0.1482
#> s3    0.1470      0.1684
#> s4    0.1759      0.1985
#> s5    0.2089      0.1581
#> s6    0.2462      0.3898
#> 
#> Lagrange multipliers (lambda; info-matrix SE, Var(lambda) = I^-1):
#>        Estimate  Std. Error  z value
#> faces    0.2068      3.1054   0.0666
#>         -0.0039      0.4124  -0.0094
#> 
#> Normalized entropy S = 0.9749   pseudo-R2 = 0.02515
#> Fano (sec 7.5): modal error pe = 0.7538 >= bound 0.588
#> Convergence code: 0 (converged) 
#> 
#> Note: the SEs above are information-matrix (curvature) quantities,
#> Var(lambda) = I^-1 (Golan 2008 sec 4.2), measuring how well the moments
#> identify lambda -- not frequentist sampling SEs (a pure inverse problem
#> is deterministic). For sampling-based inference use the
#> stochastic-moment inverse_noise().

vcov(fit2)                       # I^{-1}, the T x T dual covariance
#>           faces           
#> faces  9.643578 -1.2560837
#>       -1.256084  0.1700659
fit2$se_lambda                   # sqrt(diag(vcov)) for lambda
#>     faces           
#> 3.1054111 0.4123905
fit2$se_p                        # delta-method curvature SEs for p_hat
#>        s1        s2        s3        s4        s5        s6 
#> 0.2583316 0.1482227 0.1683946 0.1984702 0.1581484 0.3898053

5 Normalized entropy and Fano bounds

The normalized entropy summarises how much uncertainty remains after conditioning on the data. Relative to a uniform reference it is

\[\begin{equation} S(\hat{\mathbf p}) = \frac{-\sum_k \hat p_k \log(\hat p_k)}{\log(K)} \in [0, 1], \tag{5.1} \end{equation}\]

equal to \(1\) under complete ignorance (uniform \(\hat{\mathbf p}\)) and \(0\) under perfect certainty (a point mass). With a genuine prior it generalises to the prior-relative form \(S = H(\hat{\mathbf p}) / H(\mathbf{p}^0)\),

\[\begin{equation} S(\hat{\mathbf p}) = \frac{-\sum_k \hat p_k \log(\hat p_k)} {-\sum_k p_k^0 \log(p_k^0)}, \tag{5.2} \end{equation}\]

which is what inverse_ce() reports as fit$S:

c(ME = fit_me$S, CE = fit_ce$S)
#>        ME        CE 
#> 0.9005568 1.0539417

The recovered \(\hat{\mathbf p}\) is a single distribution over \(K\) states, so Golan’s Fano inequality applies directly: the modal-classification error \(p_e = 1 - \max_k \hat p_k\) is bounded below by \(S_u(\hat{\mathbf p}) - \log(2)/\log(K)\), where \(S_u = H(\hat{\mathbf p})/\log(K)\) is the uniform-reference normalized entropy. fano_bounds() returns these quantities:

fano_bounds(fit2)
#>       p_max        pe        H         S  pe_lower
#> 1 0.2461965 0.7538035 1.746698 0.9748505 0.5879977

6 From exact to noisy moments

Everything above treats the moments in (1.1) as holding exactly. When the moments are instead measured with noise — \(\mathbf{y} = \mathbf{X}\mathbf{p} + \boldsymbol\varepsilon\) — the problem becomes a generalized maximum-entropy / cross-entropy problem, and the multipliers acquire genuine sampling standard errors. That stochastic-moment sibling is inverse_noise(), which shares inverse_ce()’s formula interface and adds a noise support and signal/noise weight. As the noise support shrinks to zero, its estimates converge to the inverse_ce() solution developed here.