\(\newcommand{\B}[1]{ {\bf #1} }\) \(\newcommand{\R}[1]{ {\rm #1} }\) \(\newcommand{\W}[1]{ \; #1 \; }\)
user_predict_fit.py¶
View page sourcePredict Average Integrand Using Results of a Fit¶
Purpose¶
This examples used the fit both command to estimate the model variables. It then uses the predict fit_var command to compute the susceptible population \(S(a)\) at age \(a = 50\).
Note Table¶
north_america
/ \
united_states canada
Problem Parameters¶
The following values are used to simulate the data and set the priors for fitting the data:
iota_north_america = 1e-2
canada_effect = 0.2
united_states_effect = - canada_effect
Age and Time Values¶
The age and time values do not affect the fitting for this problem because all the functions are constant in age and time. This follows from the fact that all of the smoothings have one age and one time point.
Rate Table¶
The rate_table only specifies that the only rate variables are iota for the parent and children. In addition, it specifies the smoothings for these rates each of which has one grid point.
Variables¶
There are three model variables in this example:
iota_north_america |
|
canada_effect |
|
united_states_effect |
|
Integrand Table¶
The integrand_table for this example includes Sincidence and susceptible .
Data Table¶
There are three measurements of Sincidence in the data_table , one for each node. No noise is added to the measurements, and the priors on iota are uniform, so the fit should correspond to the model values used to simulate the data.
Avgint Table¶
There are three predictions of the susceptible population at age 50 specified in the avgint_table , one for each node.
Source Code¶
# begin problem parameters
iota_north_america = 1e-2
canada_effect = 0.2
united_states_effect = - canada_effect
# end problem parameters
# ---------------------------------------------------------------------------
import sys
import os
import copy
import math
test_program = 'example/user/predict_fit.py'
check_program = sys.argv[0].replace('\\', '/')
if check_program != test_program or len(sys.argv) != 1 :
usage = 'python3 ' + test_program + '\n'
usage += 'where python3 is the python 3 program on your system\n'
usage += 'and working directory is the dismod_at distribution directory\n'
sys.exit(usage)
print(test_program)
#
# import dismod_at
local_dir = os.getcwd() + '/python'
if( os.path.isdir( local_dir + '/dismod_at' ) ) :
sys.path.insert(0, local_dir)
import dismod_at
#
# change into the build/example/user directory
if not os.path.exists('build/example/user') :
os.makedirs('build/example/user')
os.chdir('build/example/user')
# ---------------------------------------------------------------------------
def example_db (file_name) :
def fun_iota_child(a, t) :
return ('prior_iota_child', None, None)
def fun_iota_parent(a, t) :
return ('prior_iota_parent', None, None)
# ----------------------------------------------------------------------
# age table
age_list = [ 0.0, 100.0 ]
#
# time table
time_list = [ 1995.0, 2015.0 ]
node_table = [
{ 'name':'north_america', 'parent':'' },
{ 'name':'united_states', 'parent':'north_america' },
{ 'name':'canada', 'parent':'north_america' }
]
#
# weight table:
weight_table = list()
# integrand table
integrand_table = [
{ 'name':'Sincidence' },
{ 'name':'susceptible' }
]
#
# covariate table: no covariates
covariate_table = list()
#
# mulcov table
mulcov_table = list()
#
# nslist_dict:
nslist_dict = dict()
# ----------------------------------------------------------------------
# prior_table
prior_table = [
{ # prior_iota_parent
'name': 'prior_iota_parent',
'density': 'uniform',
'lower': iota_north_america / 100.0,
'upper': iota_north_america * 100.0,
'mean': iota_north_america * 3.0,
},{ # prior_iota_child
'name': 'prior_iota_child',
'density': 'uniform',
'mean': 0.0,
}
]
# ----------------------------------------------------------------------
# smooth table
smooth_table = [
{ # smooth_iota_parent
'name': 'smooth_iota_parent',
'age_id': [ 0 ],
'time_id': [ 0 ],
'fun': fun_iota_parent
}, { # smooth_iota_child
'name': 'smooth_iota_child',
'age_id': [ 0 ],
'time_id': [ 0 ],
'fun': fun_iota_child
}
]
# ----------------------------------------------------------------------
# rate table
rate_table = [
{
'name': 'iota',
'parent_smooth': 'smooth_iota_parent',
'child_smooth': 'smooth_iota_child',
}
]
# --------------------------------------------------------------------
# option_table
option_table = [
{ 'name':'parent_node_name', 'value':'north_america' },
{ 'name':'ode_step_size', 'value':'1.0' },
{ 'name':'rate_case', 'value':'iota_pos_rho_zero' },
{ 'name':'quasi_fixed', 'value':'true' },
{ 'name':'max_num_iter_fixed', 'value':'30' },
{ 'name':'print_level_fixed', 'value':'0' },
{ 'name':'tolerance_fixed', 'value':'1e-10' },
]
# ----------------------------------------------------------------------
# data table:
data_table = list()
row = {
'density': 'log_gaussian',
'eta': '0.0',
'weight': '',
'hold_out': False,
'time_lower': 2000.0,
'time_upper': 2000.0,
'age_lower': 0.0,
'age_upper': 100.0,
'integrand': 'Sincidence',
'subgroup': 'world',
}
# north_america
row['node'] = 'north_america';
row['meas_value'] = iota_north_america
row['meas_std'] = row['meas_value'] / 10.0
data_table.append( copy.copy(row) )
# canada
row['node'] = 'canada';
row['meas_value'] = math.exp(canada_effect) * iota_north_america
row['meas_std'] = row['meas_value'] / 10.0
data_table.append( copy.copy(row) )
# united_states
row['node'] = 'united_states';
row['meas_value'] = math.exp(united_states_effect) * iota_north_america
row['meas_std'] = row['meas_value'] / 10.0
data_table.append( copy.copy(row) )
# ----------------------------------------------------------------------
# avgint table:
avgint_table = list()
# values that are the same for all data rows
row = {
'integrand': 'susceptible',
'weight': '',
'time_lower': 2000.0,
'time_upper': 2000.0,
'age_lower': 50.0,
'age_upper': 50.0,
'subgroup': 'world'
}
row['node'] = 'north_america'
avgint_table.append( copy.copy(row) )
row['node'] = 'canada'
avgint_table.append( copy.copy(row) )
row['node'] = 'united_states'
avgint_table.append( copy.copy(row) )
# ----------------------------------------------------------------------
# subgroup_table
subgroup_table = [ { 'subgroup':'world', 'group':'world' } ]
# ----------------------------------------------------------------------
# create database
dismod_at.create_database(
file_name,
age_list,
time_list,
integrand_table,
node_table,
subgroup_table,
weight_table,
covariate_table,
avgint_table,
data_table,
prior_table,
smooth_table,
nslist_dict,
rate_table,
mulcov_table,
option_table
)
# ===========================================================================
file_name = 'example.db'
example_db(file_name)
#
program = '../../devel/dismod_at'
dismod_at.system_command_prc([ program, file_name, 'init' ])
dismod_at.system_command_prc([ program, file_name, 'fit', 'both' ])
dismod_at.system_command_prc([ program, file_name, 'predict', 'fit_var' ])
# -----------------------------------------------------------------------
# connect to database
connection = dismod_at.create_connection(
file_name, new = False, readonly = True
)
predict_table = dismod_at.get_table_dict(connection, 'predict')
avgint_table = dismod_at.get_table_dict(connection, 'avgint')
node_table = dismod_at.get_table_dict(connection, 'node')
connection.close()
#
# check that all the avgint_table values were predicted (no subsetting)
assert len(predict_table) == 3
#
# S(a) = exp( - iota * a )
iota_canada = math.exp(canada_effect) * iota_north_america
iota_united_states = math.exp(united_states_effect) * iota_north_america
S_north_america = math.exp( - iota_north_america * 50.0 )
S_canada = math.exp( - iota_canada * 50.0 )
S_united_states = math.exp( - iota_united_states * 50.0 )
truth = {
'north_america' : S_north_america,
'canada' : S_canada,
'united_states' : S_united_states
}
for i in range(3) :
avgint_id = predict_table[i]['avgint_id']
node_id = avgint_table[avgint_id]['node_id']
node = node_table[node_id]['node_name']
check = truth[node]
value = predict_table[i]['avg_integrand']
assert( abs( value / check - 1.0 ) ) < 1e-6
# -----------------------------------------------------------------------------
print('predict_fit.py: OK')