sample_conditional¶

Sample Posterior for Fixed Effects Using Conditional Covariance¶

Syntax¶

mixed_object . sample_conditional (
      sample ,
      information_info ,
      solution ,
      fixed_lower ,
      fixed_upper ,
      random_opt
)

Replaced¶

This routine uses the Conditional Covariance to sample the fixed effects. This required inverting the information_mat . It has been replaced by using the sample_fixed because it is faster.

Prototype¶

void cppad_mixed::sample_conditional(
   CppAD::vector<double>&                 sample               ,
   const CppAD::mixed::sparse_mat_info&   information_info     ,
   const CppAD::mixed::fixed_solution&    solution             ,
   const CppAD::vector<double>&           fixed_lower          ,
   const CppAD::vector<double>&           fixed_upper          ,
   const CppAD::vector<double>&           random_opt           )

Purpose¶

This routine draw samples from the asymptotic posterior distribution for the optimal fixed effects (given the model and the data).

manage_gsl_rng¶

It is assumed that get_gsl_rng will return a pointer to a GSL random number generator.

mixed_object¶

We use mixed_object to denote an object of a class that is derived from the cppad_mixed base class.

sample¶

The size sample . size () is a multiple of n_fixed . The input value of its elements does not matter. We define

n_sample = sample_size / n_fixed

Upon return, for i = 0 , …, n_sample -1 , j = 0 , …, n_fixed -1 ,

sample [ i * n_fixed + j ]

is the j-th component of the i-th sample of the optimal fixed effects \(\hat{\theta}\). These samples are independent for different \(i\), and for fixed \(i\), they have the Conditional Covariance \(D\).

information_info¶

This is a sparse matrix representation for the lower triangle of the observed information matrix corresponding to solution ; i.e., the matrix returned by

information_info = mixed_object . information_mat (
      solution , random_options , random_lower , random_upper , random_in
)

solution¶

is the solution for a the call to optimize_fixed corresponding to information_info .

fixed_lower¶

is the same as fixed_lower in the call to optimize_fixed that corresponds to solution .

fixed_upper¶

is the same as fixed_upper in the call to optimize_fixed that corresponds to solution .

random_opt¶

is the optimal random effects corresponding to the solution; i.e.

      random_opt = mixed_object . optimize_random (
            random_options ,
            solution . fixed_opt ,
            random_lower ,
            random_upper ,
            random_in
      )

random_options , random_lower , random_upper , and random_in , are the same as in the call to optimize_fixed that corresponds to solution .

Theory¶

Notation¶

Given two random vectors \(u\) and \(v\), we use the notation \(\B{C}( u , v )\) for the corresponding covariance matrix; i.e.,

\[\B{C}( u , v ) = \B{E} \left( [ u - \B{E} (u) ] [ v - \B{E} (v) ]^\R{T} \right)\]

Fixed Effects Subset¶

We use \(\alpha\) for the vector of fixed effects that do not have their upper or lower bound active (or equal); i.e., if j is such that

solution . fixed_lag [ j ] == 0.0 && fixed_lower [ j ] < fixed_upper [ j ]

then \(\theta_j\) is one of the components in \(\alpha\). Note that each value of \(\alpha\) has a corresponding value for \(\theta\) where the active bounds are used for the components not in \(\alpha\).

Unconstrained Subset Covariance¶

Note that the bound constraints do not apply to the subset of fixed effects represented by \(\alpha\). We use \(\tilde{L} ( \alpha )\) to denote the fixed effects objective as a function of \(\alpha\) and where the absolute values terms in fix_likelihood are excluded. We use \(\tilde{\alpha}\) for the unconstrained optimal estimate of the subset of fixed effects and approximate its auto-covariance by

\[\B{C} ( \tilde{\alpha} , \tilde{\alpha} ) = H^{-1}\]

Here \(H\) is the Hessian corresponding to information_info . Note that information_info is the observed information matrix corresponding to all the fixed effects \(\theta\).

Constraint Equations¶

Let \(n\) be the number of fixed effects in \(\alpha\), \(m\) the number of active constraints (not counting bounds), and the equations \(e( \alpha ) = b\) be those active constraints. Here \(e : \B{R}^n \rightarrow \B{R}^m\) and \(b \in \B{R}^m\) and the inequality constraints have been converted to equalities at the active bounds (excluding the bounds on the fixed effects). Define the random variable the approximation for \(e( \alpha )\) by

\[\tilde{e} ( \alpha ) = e \left( \hat{\alpha} \right) + e^{(1)} \left( \hat{\alpha} \right) \left( \alpha - \hat{\alpha} \right)\]

where \(\hat{\alpha}\) is the subset of the optional estimate for the fixed effects solution . fixed_opt .

Conditional Covariance¶

We approximate the distribution for \(\tilde{\alpha}\) normal, and the distribution for \(\hat{\alpha}\) as the conditional distribution of \(\tilde{\alpha}\) given the value of \(\tilde{e} ( \tilde{\alpha} )\); i.e.,

\[\B{C} \left( \hat{\alpha} \W{,} \hat{\alpha} \right) = \B{C} \left( \tilde{\alpha} \W{,} \tilde{\alpha} \right) - \B{C} \left( \tilde{\alpha} \W{,} \tilde{e} \right) \B{C} \left( \tilde{e} \W{,} \tilde{e} \right)^{-1} \B{C} \left( \tilde{e} \W{,} \tilde{\alpha} \right)\]

Using the notation \(D = \B{C} \left( \hat{\alpha} \W{,} \hat{\alpha} \right)\), \(C = \B{C} \left( \tilde{\alpha} \W{,} \tilde{\alpha} \right)\), \(E = e^{(1)} \left( \hat{\alpha} \right)\), we have

\[D = C - C E^\R{T} \left( E C E^\R{T} \right)^{-1} E C\]

Example¶

The file sample_fixed.cpp is an example and test of sample_conditional was used before it was Replaced .