------------------------------------------------
lines 5-89 of file: example/user/sum_residual.py
------------------------------------------------

# {xrst_begin user_sum_residual.py}
# {xrst_comment_ch #}
#
# Sum of Residuals at Optimal Estimate
# ####################################
#
# Problem
# *******
# For this case the only data is we a sequence of  positive measurements of
# :ref:`avg_integrand@Integrand, I_i(a,t)@Sincidence`
# which we denote by :math:`y_i1` for :math:`i = 0 , \ldots , N-1`.
# We model :math:`y_i` as independent and Gaussian with mean equal
# to the true value of iota :math:`\iota`,
# and standard deviation :math:`\delta_i`.
# The negative log likelihood, up to a constant w.r.t :math:`\iota`, is
#
# .. math::
#
#   f( \iota ) =
#   \frac{1}{2} \sum_{i=0}^{N-1} \left( \frac{y_i - \iota}{\delta_i} \right)^2
#
# Optimal Solution
# ****************
# The optimal estimator for :math:`\iota` satisfies the equation
#
# .. math::
#
#   0 = f'( \hat{\iota} ) =
#       - \sum_{i=1}^{N-1} \frac{y_i - \hat{\iota} }{\delta_i^2}
#
# Weighted Residuals
# ******************
# We use the notation :math:`r_i` for the weighted residuals
#
# .. math::
#
#   r_i = \frac{y_i - \hat{\iota}}{\delta_i}
#
# If :math:`\hat{\iota}` were the true value for :math:`\iota`,
# the weighted residuals would be mean zero and variance one.
# But :math:`\hat{\iota}` is instead the maximum likelihood estimator and
#
# .. math::
#
#   0 = \sum_{i=1}^{N-1} \frac{y_i - \hat{\iota} }{\delta_i^2}
#
# .. math::
#
#   0 = \sum_{i=1}^{N-1} \frac{r_i}{\delta_i}
#
# Note that if :math:`\delta_i` were the same for all :math:`i`,
# the sum of the weighted residuals :math:`\sum_i r_i` would be zero.
#
# CV Standard Deviations
# **********************
# We consider the case were we a coefficient of variation :math:`\lambda`
# is used to model the measurement noise; :math:`\delta_i = \lambda y_i`.
# In this case
#
# .. math::
#
#   0 = \sum_{i=1}^{N-1} \frac{r_i}{y_i}
#
# Weighted Average of Weighted Residuals
# **************************************
# We define the weight :math:`w_i` by
#
# .. math::
#
#   w_i = \delta_i^{-1} / \sum_{i=0}^{N-1} \delta_i^{-1}
#
# The corresponding weighted average of the weighted residuals is zero; i.e,
#
# .. math::
#
#   0 = \sum_{i=1}^{N-1} w_i r_i
#
# Source Code
# ***********
# {xrst_literal
#       BEGIN PYTHON
#       END PYTHON
# }
#
# {xrst_end user_sum_residual.py}