-------------------------------------------------
lines 5-191 of file: example/user/one_function.py
-------------------------------------------------

# {xrst_begin user_one_function.py}
# {xrst_spell
#       ds
# }
# {xrst_comment_ch #}
#
# Fitting One Function of Two Variables
# #####################################
#
# Purpose
# *******
# This example shows how to fit one function of two variables given
# direct measurements of the function.
# For our example, we are give measurements of prevalence,
# as a function of age and time, and fit a function to the measurements.
#
# rho
# ***
# We use *rho*
# as a surrogate for prevalence because our model for prevalence
# does not use *rho* (we could have used any other rate
# and its corresponding integrand).
#
# Random Effects
# **************
# To keep this example simple, there is only one node ``world``
# in the :ref:`node_table-name` , and the
# :ref:`mulcov_table@subgroup_smooth_id`
# is null in the mulcov_table.
# Hence, there are not random effects in this example.
#
# Covariates
# **********
# Income is the only covariate for this example and it has been normalized
# to be between zero and one.  The income reference value corresponds to
# no effect. The income multiplier is the true value used to simulate the
# data.
# {xrst_spell_off}
# {xrst_code py}
income_reference   =   0.5
income_multiplier  = - 0.2
income_mulcov_type = "meas_value"
# {xrst_code}
# {xrst_spell_on}
# The *income_mulcov_type* can be either ``"meas_value"`` or
# ``"rate_value"`` .
# It should not make a difference in the results which
# :ref:`mulcov_table@mulcov_type` is used because
# there is only one rate (one function) being fit, and the
# :ref:`integrand_table@integrand_name@ODE` is not needed to compute
# the integrand.
# You can test this by changing *income_mulcov_type* to
# ``"rate_value"``
# and uses the same :ref:`user_one_function.py@random_seed` .
#
# Simulated Data
# **************
#
# Data Density
# ============
# A direct measurement of *rho*
# (which we are using to represent prevalence) corresponds to the
# :ref:`remission integrand<avg_integrand@Integrand, I_i(a,t)@remission>` .
# We use a :ref:`density_table@density_name@log_gaussian`
# model with the following parameters:
# {xrst_spell_off}
# {xrst_code py}
prevalence_delta  = 0.1
prevalence_eta    = 1e-6
# {xrst_code}
# {xrst_spell_on}
#
# True Prevalence
# ===============
# For simulating our prevalence data we consider the case where
# all the rates are constant and *rho* , *chi* are zero; i.e.,
# the differential equation is
#
# .. math::
#   :nowrap:
#
#   \begin{eqnarray}
#   S'(a) = -( \iota + \omega ) S(a)
#   \\
#   C'(a) = + \iota S(a) - \omega C(a)
#   \end{eqnarray}
#
# where :math:`S(0) = 1` and :math:`C(0) = 0`.
# It follows that
#
# .. math::
#
#   S(a) = \exp[ - ( \iota + \omega ) a ]
#
# .. math::
#
#   C(a)
#   = \int_0^a \exp[ \omega (s - a) ] \iota S(s) ds
#
# You can check the formula for :math:`C(a)` as follows:
# differentiating w.r.t. :math:`a` inside the integral yields
# :math:`- \omega C(a)` and differentiating the upper limit
# w.r.t. :math:`a` yields :math:`\iota S(a)`.
# It follows that
#
# .. math::
#
#   C(a) = \iota \int_0^a \exp[ \omega (s - a) - ( \iota + \omega) s ] ds
#
# .. math::
#
#   C(a) = \iota \exp( - \omega a) \int_0^a \exp( - \iota s ) ds
#
# .. math::
#
#   C(a) = \exp( - \omega a) - \exp[ - ( \iota + \omega) a ]
#
# {xrst_spell_off}
# {xrst_code py}
def Prevalence(age) :
    # rho and chi are zero for this simulation of prevalence data
    iota      = 1e-4 # incidence used to simulate prevalence data
    omega     = 1e-2 # other cause mortality used to simulate prevalence data
    exp_omega = exp(-omega * age)
    exp_sum   = exp(-(iota + omega) * age)
    S         = exp_sum
    C         = exp_omega - exp_sum
    return C / (S + C)
# {xrst_code}
# {xrst_spell_on}
#
# Computing Delta
# ***************
# Once we have simulated a measurement value,
# we can solve for :math:`\Delta` are follows; see
# :ref:`delta<data_like@Transformed Standard Deviation, delta_i(theta)>` :
#
# .. math::
#
#   \delta = \log( y + \eta + \Delta ) - \log(y + \eta )
#
# .. math::
#
#   \exp ( \delta ) = \frac{ y + \eta + \Delta }{ y + \eta }
#
# .. math::
#
#   \exp ( \delta ) ( y + \eta ) = y + \eta + \Delta
#
# .. math::
#
#   \Delta = [ \exp ( \delta ) - 1 ] ( y + \eta )
#
# For this case there are no measurement noise covariates so
# :math:`\delta` is the standard deviation for the simulated data.
# Furthermore, the
# :ref:`option_table@minimum_meas_cv` is zero,
# so :math:`\Delta` is the *meas_std* .
#
# Prevalence Prior
# ****************
# Prevalence must always be between zero and one so this limits are used
# in the prior for prevalence.
# Some functions might allow for negative values and the lower limit
# for the rate would be negative in that case.
#
# random_seed
# ***********
# We use the clock to choose a seed for the random number generator.
# If an error occurs, the seed will be printed so that the error can be
# reproduced. You can also use a fixed value for the random seed to see
# how changing other parameters changes the results.
# {xrst_spell_off}
# {xrst_code py}
import time
random_seed = int( time.time() )
# {xrst_code}
# {xrst_spell_on}
#
# Source Code
# ***********
# {xrst_literal
#       BEGIN PYTHON
#       END PYTHON
# }
#
# {xrst_end user_one_function.py}