csv.predict_descend¶

Example Predicting From Self and a Descendant¶

iota_true¶

There are only two non-zero rate in this model, omega and the parent value for iota (the model for incidence).

omega_true = 0.01
iota_true  = 0.01

node_dict¶

Keys are nodes and values are corresponding parent node:

node_dict = {
    'n0' : ''   ,
    'n1' : 'n0' ,
    'n2' : 'n0' ,
}

age_grid, time_grid¶

The only rate for this example is iota and it is constant in age and time so we use a very simple age and time grid:

age_grid   = [0.0, 100.0]
time_grid  = [1980.0, 2020.0]

number_sample¶

This is the number of samples of the posterior distribution to generate for each successful fit. Note that each sample includes values for all the dismod_at model variables-name . For this example there is only one model variable iota .

number_sample = 4000

descendant_std_factor¶

This is the descendant_std_factor . This example checks that this factor is (is not) used when predicting for the node that was fit (for a descendant of the node that was fit).

descendant_std_factor = 2.0

random_seed¶

We us the current time in seconds to seed the random number generator:

import time
random_seed = int( time.time() )

n_data, std_data¶

There are n_data simulated measurements. These measurements are a Gaussian with mean iota_true and standard deviation std_data :

n_data = 1
std_data = iota_true * 0.2

input_file¶

Input CSV files that are placed in the fit directory:

input_file = dict()

predict_integrand.csv¶

All the measurements and predictions for this example are for Sincidence; i.e.. a direct measurement of iota:

input_file['predict_integrand.csv']  = 'integrand_name\n'
input_file['predict_integrand.csv'] += 'Sincidence\n'

There are measurements for n0, but not for n1 or n2. Thus node n0 will predict for itself and for n2 ( n1 is not in fit_goal.csv . )

option_fit.csv¶

input_file['option_fit.csv']  =  \
"""name,value
refit_split,false
"""
input_file['option_fit.csv'] += f'number_sample,{number_sample}\n'
input_file['option_fit.csv'] += f'random_seed,{random_seed}\n'

option_predict.csv¶

input_file['option_predict.csv']  =  'name,value\n'
input_file['option_predict.csv']  +=  \
    f'descendant_std_factor,{descendant_std_factor}\n'

fit_goal.csv¶

There will be no fit or predictions for n1:

input_file['fit_goal.csv']  = 'node_name\n'
input_file['fit_goal.csv'] += 'n2'

node.csv¶

input_file['node.csv'] = \
'node_name,parent_name\n'
for node_name in node_dict :
    parent_name = node_dict[node_name]
    input_file['node.csv'] += f'{node_name},{parent_name}\n'

covariate.csv¶

input_file['covariate.csv'] = 'node_name,sex,age,time,omega\n'
for node_name in [ 'n0', 'n1', 'n2' ] :
    for sex in [ 'female', 'male' ] :
        for age in age_grid :
            for time in time_grid :
                row   = f'{node_name},{sex},{age},{time},{omega_true}\n'
                input_file['covariate.csv'] += row

prior.csv¶

uniform_eps_0.1 is uniform on the interval [eps,.1] where eps = 1e-6, the mean, which is only used for initialization, is 0.03.

input_file['prior.csv'] = \
"""name,density,mean,std,eta,lower,upper
uniform_eps_0.1,uniform,0.01,,,1e-6,0.1
"""

parent_rate.csv¶

The parent rate for iota is has a uniform_eps_0.1 prior and is constant w.r.t. age and time:

input_file['parent_rate.csv'] = \
"""rate_name,age,time,value_prior,dage_prior,dtime_prior,const_value
iota,0.0,0.0,uniform_eps_0.1,,,
"""

child_rate.csv¶

All the child rates (rate random effects) are zero

input_file['child_rate.csv'] = 'rate_name,value_prior\n'

mulcov.csv¶

There are no covariate multipliers

input_file['mulcov.csv'] = 'covariate,type,effected,value_prior,const_value\n'

Hessian of Likelihood¶

For this example, iota is the only model variable (that does not have equal lower and upper bounds). Let \(m\) be n_data , \(sigma\) be std_data, and \(y_i\) be the measurement values. The negative log likelihood for the n0 fit is:

\[L( \iota , y ) = \frac{1}{2} \sum_{i=0}^{m} \log \left( 2 \pi \sigma^2 \right) + \left( \frac{y_i - \iota }{\sigma} \right)^2\]

The Hessian with respect to iota ( \(\iota\) ) is

\[\frac{ \partial^2 }{\partial \iota^2} L = \frac{m}{ \sigma^2 }\]

Rest of Source Code¶

Below is the source code, except for the settings above.

#
# os, sys, copy
import os
import sys
import copy
#
# random
import random
random.seed(random_seed)
#
# at_cascade
current_directory = os.getcwd()
if os.path.isfile( current_directory + '/at_cascade/__init__.py' ) :
    sys.path.insert(0, current_directory)
import at_cascade
#
# fit
def create_input_files(fit_dir) :
    #
    # csv files in input_file
    for name in input_file :
        file_name = f'{fit_dir}/{name}'
        file_ptr  = open(file_name, 'w')
        file_ptr.write( input_file[name] )
        file_ptr.close()
    #
    # data_table
    data_table = list()
    data_id    = -1
    age        = 50.0
    time       = 2000.0
    row = {
        'node_name'      : 'n0',
        'integrand_name' : 'Sincidence',
        'hold_out'       : 0,
        'density_name'   : 'gaussian',
        'age_lower'      : age,
        'age_upper'      : age,
        'time_lower'     : time,
        'time_upper'     : time,
        'eta'            : None,
        'nu'             : None,
        'meas_std'       : std_data,
    }
    for i in range(n_data) :
        if i % 2 == 0 :
            sex = 'female'
        else :
            sex = 'male'
        data_id          += 1
        row['data_id']    = data_id
        row['sex']        = sex
        row['meas_value'] = random.normalvariate( iota_true, std_data )
        data_table.append( copy.copy( row ) )
    #
    # data_in.csv
    at_cascade.csv.write_table(
        file_name = f'{fit_dir}/data_in.csv' ,
        table     = data_table ,
    )
#
# main
def main() :
    #
    # fit_dir
    fit_dir = 'build/example/csv/fit'
    at_cascade.empty_directory(fit_dir)
    #
    # create_input_files
    create_input_files(fit_dir)
    #
    # at_cascade.csv.fit
    at_cascade.csv.fit(fit_dir)
    #
    # at_cascade.csv.predict
    at_cascade.csv.predict(fit_dir)
    #
    # fit
    file_name = f'{fit_dir}/fit_predict.csv'
    fit_predict_table = at_cascade.csv.read_table(file_name);
    fit               = dict()
    for row in fit_predict_table :
        avgint_id = int( row['avgint_id'] )
        node_name = row['node_name']
        sex       = row['sex']
        key       = (avgint_id, node_name, sex)
        assert key not in fit
        fit[key]  = float( row['avg_integrand']  )
    #
    # sam
    file_name = f'{fit_dir}/sam_predict.csv'
    sam_predict_table = at_cascade.csv.read_table(file_name);
    sam               = dict()
    for row in sam_predict_table :
        avgint_id = int( row['avgint_id'] )
        node_name = row['node_name']
        sex       = row['sex']
        key       = (avgint_id, node_name, sex)
        if key not in sam :
            sam[key] = list()
        sam[key].append(  float( row['avg_integrand']  ) )
    #
    # hessian_like
    hessian_like = n_data / ( std_data * std_data )
    #
    # check predctions for n0
    for avgint_id in range(4) :
        key   = ( avgint_id, 'n0', 'both' )
        sumsq = 0.0
        assert number_sample == len( sam[key] )
        for avg_integrand in sam[key] :
            sumsq += ( avg_integrand - fit[key] )**2
        sam_cov = sumsq / number_sample
        relerr  = 1.0 -  sam_cov * hessian_like
        assert abs( relerr ) < 0.05
    #
    # check predctions for n2
    for avgint_id in range(4) :
        for sex in [ 'female', 'male' ] :
            key   = ( avgint_id, 'n2', sex )
            sumsq = 0.0
            assert number_sample == len( sam[key] )
            for avg_integrand in sam[key] :
                sumsq  += ( avg_integrand - fit[key] )**2
            sam_cov    =  sumsq / number_sample
            cov_factor = descendant_std_factor * descendant_std_factor
            relerr     = 1.0 -  (sam_cov / cov_factor) * hessian_like
            assert abs( relerr ) < 0.05
    #
    #
    print('predict_descend.py: OK')
#
if __name__ == '__main__' :
    main()
#