-----------------------------------------------
lines 5-120 of file: example/user/hes_random.py
-----------------------------------------------

# {xrst_begin user_hes_random.py}
# {xrst_spell
#       cc
# }
# {xrst_comment_ch #}
#
# Check the Hessian of the Random Effects Objective
# #################################################
#
# Purpose
# *******
# This example checks the computation of the Hessian of the
# random effects objective used by the
# :ref:`sample_command@asymptotic` sampling method.
#
# Reference
# *********
# See the
# `theory <https://bradbell.github.io/cppad_mixed/doc/theory.htm>`_
# section of the ``cppad_mixed`` documentation.
#
# Notation
# ********
#
# .. csv-table::
#       :widths: auto
#
#       :math:`\theta`,model incidence for the parent region
#       :math:`u_0`,incidence random effect for first child
#       :math:`u_1`,incidence random effect for second child
#       :math:`y_0`,measured incidence for the first child
#       :math:`y_1`,measured incidence for the second child
#       :math:`s`,standard deviation for data and random effects
#
# The only fixed effect in this model is :math:`\theta`
# (sometimes written *theta* ) the incidence level for the world.
# The random effects are :math:`u_0` and :math:`u_1`.
#
# Problem Parameters
# ******************
# {xrst_spell_off}
# {xrst_code py}
import csv
import math
import time
theta_true     = 0.5
u_true         = [ 0.5, -0.5 ]
y_0            = theta_true * math.exp( u_true[0] )
y_1            = theta_true * math.exp( u_true[1] )
standard_dev   = theta_true / 10.0 # the standard deviation s
random_seed    = int( time.time() )
number_sample  = 4000
# {xrst_code}
# {xrst_spell_on}
#
# Random Likelihood
# *****************
# The negative log-likelihood for this example, ignoring the constant
# of integration, is
#
# .. math::
#
#   f( \theta, u )
#   = \frac{1}{2 s^2} \left(
#       [ y_0 - \theta \exp( u_0 ) ]^2 +
#       [ y_1 - \theta \exp( u_1 ) ]^2 +
#       u_0^2 +
#       u_1^2
#   \right)
#
# Gradient w.r.t. Random Effects
# ******************************
#
# .. math::
#
#   f_u ( \theta, u )
#   =
#   \frac{1}{s^2} \left(
#   \begin{array}{c}
#       \theta^2 \exp( 2 u_0 ) - y_0 \theta \exp( u_0 )  + u_0
#       \\
#       \theta^2 \exp( 2 u_1 ) - y_1 \theta \exp( u_1 )  + u_1
#   \end{array}
#   \right)
#
# Hessian w.r.t. Random Effects
# *****************************
#
# .. math::
#
#   f_{u,u} ( \theta, u )
#   =
#   \frac{1}{s^2} \left(
#   \begin{array}{cc}
#       2 \theta^2 \exp( 2 u_0 ) - y_0 \theta \exp( u_0 ) + 1   & 0
#       \\
#       0   & 2 \theta^2 \exp( 2 u_1 ) - y_1 \theta \exp( u_1 ) + 1
#   \end{array}
#   \right)
#
# Asymptotic Statistics
# *********************
# The asymptotic posterior distribution for random effects is normal
# with mean :math:`\hat{u} ( \theta )`
# and variance :math:`f_{u,u} [ \theta , \hat{u} ( \theta ) ]^{-1}`
# where :math:`\hat{u} ( \theta )` minimizes the random effects objective
# :math:`f( \theta , \cdot )`.
#
# Source Code
# ***********
# {xrst_literal
#       BEGIN PYTHON
#       END PYTHON
# }
#
# {xrst_end user_hes_random.py}