sample_random.cpp

View page source

Sample From Fixed Effects Posterior: Example and Test

# include <cppad/cppad.hpp>
# include <cppad/mixed/cppad_mixed.hpp>
# include <cppad/mixed/manage_gsl_rng.hpp>
# include <Eigen/Dense>

namespace {
   using CppAD::vector;
   using CppAD::log;
   using CppAD::AD;
   //
   using CppAD::mixed::d_sparse_rcv;
   using CppAD::mixed::d_vector;

   class mixed_derived : public cppad_mixed {
   private:
# ifndef NDEBUG
      const size_t          n_fixed_;
# endif
      const size_t          n_random_;
      const d_vector&       y_;
   // ----------------------------------------------------------------------
   public:
      // constructor
      mixed_derived(
         size_t                 n_fixed       ,
         size_t                 n_random      ,
         bool                   quasi_fixed   ,
         bool                   bool_sparsity ,
         const d_vector&       y              ) :
         cppad_mixed(
            n_fixed, n_random, quasi_fixed, bool_sparsity
         )                     ,
# ifndef NDEBUG
         n_fixed_(n_fixed)     ,
# endif
         n_random_(n_random)   ,
         y_(y)
      {  assert( n_fixed == 2);
         assert( y_.size() == n_random_ );
      }
   // ----------------------------------------------------------------------
      // implementation of ran_likelihood
      a1_vector ran_likelihood(
         const a1_vector&         theta  ,
         const a1_vector&         u      ) override
      {
         assert( theta.size() == n_fixed_ );
         assert( u.size() == y_.size() );
         a1_vector vec(1);

         // initialize part of log-density that is always smooth
         vec[0] = 0.0;

         // sqrt_2pi = CppAD::sqrt(8.0 * CppAD::atan(1.0) );

         for(size_t i = 0; i < n_random_; i++)
         {  a1_double mu     = u[i] + theta[0];
            a1_double sigma  = theta[1];
            a1_double res    = (y_[i] - mu) / sigma;

            // p(y_i | u, theta)
            vec[0] += log(sigma) + res * res / 2.0;
            // following term does not depend on fixed or random effects
            // vec[0] += log(sqrt_2pi);

            // p(u_i | theta)
            a1_double sq = u[i] * u[i];
            if( i == 0 )
               sq = (u[0] + u[1]) * (u[0] + u[1]);
            vec[0] += sq / 2.0;
            // following term does not depend on fixed or random effects
            // vec[0] += log(sqrt_2pi);
         }
         return vec;
      }
   };
}

bool sample_random_xam(void)
{
   bool   ok = true;
   double inf = std::numeric_limits<double>::infinity();
   //
   // initialize gsl random number generator
   size_t random_seed = CppAD::mixed::new_gsl_rng(0);
   //
   size_t n_data   = 10;
   size_t n_fixed  = 2;
   size_t n_random = n_data;
   //
   d_vector data(n_data), random_in(n_random);
   for(size_t i = 0; i < n_data; i++)
   {  data[i]       = double(i + 1);
      random_in[i]    = 0.0;
   }
   d_vector fixed_vec(n_fixed);
   for(size_t i = 0; i < n_fixed; i++)
      fixed_vec[i] = double(i + 1);

   // object that is derived from cppad_mixed
   bool quasi_fixed   = true;
   bool bool_sparsity = true;
   mixed_derived mixed_object(
      n_fixed, n_random, quasi_fixed, bool_sparsity, data
   );
   mixed_object.initialize(fixed_vec, random_in);

   std::string random_ipopt_options =
      "Integer print_level     0\n"
      "String  sb              yes\n"
      "String  derivative_test second-order\n"
      "Numeric tol             1e-8\n"
   ;
   d_vector random_lower(n_random), random_upper(n_random);
   for(size_t i = 0; i < n_random; i++)
   {  random_lower[i] = -inf;
      random_upper[i] = +inf;
   }
   //
   // compute the optimal random effects
   d_vector random_opt = mixed_object.optimize_random(
      random_ipopt_options, fixed_vec, random_lower, random_upper, random_in
   );
   //
   // sample from the posterior for random effects given fixed effects
   // and compute the  sample covariance matrix
   double cov_factor = 4.0;
   size_t n_sample   = 20000;
   d_vector sample(n_sample * n_random);
   std::string error_msg = mixed_object.sample_random(
      sample,
      random_ipopt_options,
      fixed_vec,
      random_lower,
      random_upper,
      random_in,
      cov_factor
   );
   ok &= error_msg == "";
   d_vector sample_cov(n_random * n_random);
   for(size_t i = 0; i < n_random; i++)
      for(size_t j = 0; j < n_random; j++)
         sample_cov[i * n_random + j] = 0;
   for(size_t i_sample = 0; i_sample < n_sample; i_sample++)
   {  d_vector diff(n_random);
      for(size_t j = 0; j < n_random; j++)
         diff[j] = sample[i_sample * n_random + j] - random_opt[j];
      for(size_t i = 0; i < n_random; i++)
         for(size_t j = 0; j < n_random; j++)
            sample_cov[i * n_random + j] +=
               diff[i] * diff[j] / double(n_sample);
   }
   //
   // For i > 1, the terms involving u[i] is
   //    0.5 * ( (y[i] - theta[0] - u[i])^2 / theta[1]^2 + u[i]^2  )
   // The terms involving u[0] and u[1] are
   //    0.5 * ( (y[1] - theta[0] - u[1])^2 / theta[1]^2 + u[1]^2  )
   //  + 0.5 * ( (y[0] - theta[0] - u[0])^2 / theta[1]^2 + (u[0]+u[1])^2  )
   double diag = 1.0 / (fixed_vec[1] * fixed_vec[1]) + 1.0;
   typedef Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > double_mat;
   double_mat info_mat = double_mat::Zero(n_random, n_random);
   for(size_t i = 0; i < n_random; i++)
      info_mat(i, i) = diag;
   info_mat(0, 1)  = 1.0;
   info_mat(1, 0)  = 1.0;
   info_mat(1, 1) += 1.0;
   //
   double_mat cov_mat = cov_factor * info_mat.inverse();
   for(size_t i = 0; i < n_random; i++)
   {  for(size_t j = 0; j < n_random; j++)
      {  double value = sample_cov[i * n_random + j];
         double check = cov_mat(i, j);
         double scale = std::sqrt( cov_mat(i,i) * cov_mat(j,j) );
         double relerr = (value - check) / scale;
         ok           &= std::fabs(relerr) < 0.05;
      }
   }
   //
   if( ! ok )
      std::cout << "\nrandom_seed = " << random_seed << "\n";
   //
   CppAD::mixed::free_gsl_rng();
   return ok;
}