-------------------------------------------- lines 15-322 of file: example/csv/no_data.py -------------------------------------------- {xrst_begin csv.no_data} {xrst_spell diota const dage delim dtime eps meas sincidence std } Use a No Data Example to Understand Priors ########################################## csv_file ******** This dictionary is used to hold the data corresponding to the csv files for this example: {xrst_code py}""" csv_file = dict() """{xrst_code} node.csv ******** For this example the root node, n0, has no children. {xrst_code py}""" csv_file['node.csv'] = \ '''node_name,parent_name n0, ''' """{xrst_code} option_fit.csv ************** This example uses the default value for the options that are not listed below. For example, refit_split is true so that we get predictions for the female / male split at node n0. (The separate female / male fits will abort because there is no data for either of these cases and they are not the root job case.) The number of samples is huge so that the sample covariance is close to the true covariance (for testing). .. csv-table:: :header: Name, Value :delim: | tolerance_fixed | this is set small, 1e-8, so we can check accuracy. random_seed | chosen using current seconds reported by python time package. number_sample | The number of samples of the posterior for each fit {xrst_code py}""" random_seed = str( int( time.time() ) ) number_sample = 10000 csv_file['option_fit.csv'] = 'name,value\n' csv_file['option_fit.csv'] += 'tolerance_fixed,1e-8\n' csv_file['option_fit.csv'] += f'random_seed,{random_seed}\n' csv_file['option_fit.csv'] += f'number_sample,{number_sample}\n' """{xrst_code} option_predict.csv ****************** This example uses the default value for all the options in option_predict.csv. (Future versions of this example will discuss the predictions.) {xrst_code py}""" csv_file['option_predict.csv'] = 'name,value\n' csv_file['option_predict.csv'] += 'db2csv,true\n' """{xrst_code} covariate.csv ************* This example has no covariates {xrst_code py}""" csv_file['covariate.csv'] = \ '''node_name,sex,age,time,omega n0,female,0,2000,0.02 n0,female,100,2000,0.02 n0,male,0,2000,0.02 n0,male,100,2000,0.02 ''' """{xrst_code} fit_goal.csv ************ This example only fits the root node n0. {xrst_code py}""" csv_file['fit_goal.csv'] = \ '''node_name n0 ''' """{xrst_code} predict_integrand.csv ********************* For this example we want to know the values of Sincidence. (Note that Sincidence is a direct measurement of iota.) {xrst_code py}""" csv_file['predict_integrand.csv'] = \ '''integrand_name Sincidence ''' """{xrst_code} iota_mean ********* This is the mean value of iota in the gauss_eps_1 prior: {xrst_code py}""" iota_mean = 0.01 """{xrst_code} iota_std ******** This is the standard deviation value of iota in the gauss_eps_1 prior: {xrst_code py}""" iota_std = 0.1 """{xrst_code} diota_std ********* This is the standard deviation of the iota age difference in gauss_0d prior: {xrst_code py}""" diota_std = 0.2 """{xrst_code} prior.csv ********* We define three priors: .. csv-table:: :widths: auto :delim: ; gauss_eps_1; a Gaussian with lower eps, upper 1, iota_mean, iota_std gauss_0_d; a Gaussian with mean 0 standard deviation diota_std {xrst_code py}""" csv_file['prior.csv'] = \ 'name,lower,upper,mean,std,density\n' \ f'gauss_eps_1,1e-6,1.0,{iota_mean},{iota_std},gaussian\n' \ f'gauss_0_d,,,0.0,{diota_std},gaussian\n' """{xrst_code} parent_rate.csv *************** The only non-zero rates are omega and iota (omega is known and specified by the covariate.csv file). The model for iota is linear w.r.t age and constant w.r.t. time. Its value prior is gauss_eps_1 and its dage prior is gauss_0_d. It does not have any dtime priors because there are no time differences between grid values. {xrst_code py}""" csv_file['parent_rate.csv'] = \ '''rate_name,age,time,value_prior,dage_prior,dtime_prior,const_value iota,0.0,2000,gauss_eps_1,gauss_0_d,, iota,100.0,2000,gauss_eps_1,gauss_0_d,, ''' """{xrst_code} child_rate.csv ************** The are no children (hence no random effects) in this example. {xrst_code py}""" csv_file['child_rate.csv'] = \ '''rate_name,value_prior ''' """{xrst_code} mulcov.csv ********** There are no covariate multipliers in this example: {xrst_code py}""" csv_file['mulcov.csv'] = \ '''covariate,type,effected,value_prior,const_value ''' """{xrst_code} data_in.csv *********** There is no data in this example {xrst_code py}""" header = 'data_id,integrand_name,node_name,sex,age_lower,age_upper,' header += 'time_lower,time_upper,meas_value,meas_std,hold_out,' header += 'density_name,eta,nu' csv_file['data_in.csv'] = header + '\n' r"""{xrst_code} Negative Log Likelihood *********************** The standard deviation in the difference prior above is one. Consider the root job; i.e., the fit for node n0 and both sexes. Two times the negative log likelihood for this example :math:`2 L(x)` is: .. math:: 2 L( x ) = - \log( 2 \pi d^2 ) + \left( \frac{ x_1 - x_0 }{ d } \right)^2 + \sum_{j=0}^1 - \log( 2 \pi \sigma^2 ) + \left( \frac{ x_j - \bar{\iota} }{ \sigma } \right)^2 where :math:`x_0` is iota at age 0, :math:`x_1` is iota at age 100, :math:`\bar{\iota}` is *iota_mean* , :math:`\sigma` is *iota_std* , :math:`d` is *diota_std* , Gradient ======== The partial of :math:`L(x)` with respect to :math:`x_j` is: .. math:: \frac{ \partial L(x)} { \partial x_j } = \frac{x_j - x_{1-j}}{ d^2 } + \frac{ x_j - \bar{\iota} }{ \sigma^2 } Hessian ======= The second partial of :math:`L(x)` with respect to :math:`x_j` is: .. math:: \frac{ \partial^2 L(x)} { \partial x_j^2 } = d^{-2} + \sigma^{-2} The cross partial with respect to :math:`x_0` and :math:`x_1` is: .. math:: \frac{ \partial^2 L(x)} { \partial x_0 \partial x_1 } = - d^{-2} The Hessian of :math:`L(x)` is: .. math:: H & = \begin{bmatrix} d^{-2} + \sigma^{-2} & -d^{-2} \\ -d^{-2} & d^{-2} + \sigma^{-2} \end{bmatrix} \\ H & = d^{-2} \begin{bmatrix} 1 + ( d / \sigma )^2 & - 1 \\ - 1 & 1 + ( d / \sigma )^2 \end{bmatrix} The determinant of :math:`d^2 H` is: .. math:: \det(d^2 H ) = 2 ( d / \sigma )^2 + ( d / \sigma )^4 The inverse of :math:`d^2 H` is: .. math:: d^{-2} H^{-1} & = \frac{1}{ 2 ( d / \sigma )^2 + ( d / \sigma )^4 } \begin{bmatrix} 1 + ( d / \sigma )^2 & 1 \\ 1 & 1 + ( d / \sigma )^2 \end{bmatrix} \\ d^{-2} H^{-1} & = \frac{ ( d / \sigma )^2 }{ 2 + ( d / \sigma )^2 } \begin{bmatrix} 1 + ( d / \sigma )^2 & 1 \\ 1 & 1 + ( d / \sigma )^2 \end{bmatrix} It follows that: .. math:: H^{-1} & = \frac{ \sigma^{-2} }{ 2 + ( d / \sigma )^2 } \begin{bmatrix} 1 + ( d / \sigma )^2 & 1 \\ 1 & 1 + ( d / \sigma )^2 \end{bmatrix} \\ H^{-1} & = \frac{ \sigma^{-2} }{ 2 ( \sigma / d )^{-2} + 1 } \begin{bmatrix} ( \sigma / d )^{-2} + 1 & ( \sigma / d )^{-2} \\ ( \sigma / d )^{-2} & ( \sigma / d )^{-2} + 1 \end{bmatrix} Covariance of sam_predict.csv ============================= The covariance of the samples in o :ref:`csv.fit@Output Files@sam_predict.csv` is equal to :math:`H^{-1}` . Note that as :math:`\sigma / d` (i.e. *iota_std* / *diota_std* ) gets small, the variance of the estimates for iota gets close to :math:`\sigma^2` (i.e., the square of *iota_std* ) . Source Code *********** {xrst_literal BEGIN_PROGRAM END_PROGRAM } {xrst_end csv.no_data}