Sample Posterior for Fixed Effects Using Conditional Covariance

Processing math: 0%

sample_conditional

$\newcommand{\R}[1]{ {\rm #1} } \newcommand{\B}[1]{ {\bf #1} } \newcommand{\W}[1]{ \; #1 \; }$ This is cppad_mixed--20220519 documentation: Here is a link to its current documentation . Sample Posterior for Fixed Effects Using Conditional Covariance
Syntax
Replaced
Prototype
Purpose
manage_gsl_rng
mixed_object
sample
information_info
solution
fixed_lower
fixed_upper
random_opt
Theory
     Notation
     Fixed Effects Subset
     Unconstrained Subset Covariance
     Constraint Equations
     Conditional Covariance
Example

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 .

Input File: src/eigen/sample_conditional.cpp