\(\newcommand{\B}[1]{ {\bf #1} }\) \(\newcommand{\R}[1]{ {\rm #1} }\) \(\newcommand{\W}[1]{ \; #1 \; }\)
user_ill_condition.py#
View page sourceAn Ill Conditioned Example Where Re-Scaling is Helpful#
Notation#
\(\bar{\chi}\) |
measured value for excess mortality at age zero |
\(\sigma\) |
std for the measurement noise |
\(\delta\) |
std for age difference for \(\chi\) |
\(\chi_0\) |
estimate for \(\chi\) at age zero |
\(\chi_1\) |
estimate for \(\chi\) at age 100 |
Objective#
For this case the negative log likelihood is, not counting terms that are constant w.r.t \(( \chi_0, \chi_1 )\),
\[f( \chi_0 , \chi_1 )
=
\frac{1}{2} \left( \frac{\chi_0 - \bar{\chi}}{\sigma} \right)^2
+
\frac{1}{2} \left( \frac{\chi_1 - \chi_0}{\delta} \right)^2\]
Derivative#
The partial of the objective w.r.t. \(\chi_0\) is
\[\partial_0 f ( \chi_0 , \chi_1 )
=
\frac{\chi_0 - \bar{\chi}}{ \sigma^2 }
-
\frac{\chi_1 - \chi_0}{\delta^2}\]
The partial of the objective w.r.t. \(\chi_1\) is
\[\partial_1 f ( \chi_0 , \chi_1 )
=
\frac{\chi_1 - \chi_0}{\delta^2}\]
Ill-Conditioning#
If \(\sigma\) is much smaller than \(\delta\), then the optimization problem will weight the measurement equation \(\chi_0 - \bar{\chi} = 0\) much more than the difference equation \(\chi_1 - \chi_0 = 0\). If \(\sigma\) is much larger than \(\delta\), then the optimization problem will weight the difference equation much more than the measurement equation.
Scaling#
Note that the initial scaling, in the example below, uses the mean values for \(\chi_0\), \(\chi_1\) which are zero.
Source Code#
# You can changed the values below and rerun this program
chi_bar = 0.01
sigma = 1e-5 * chi_bar
delta = 1e-1 * chi_bar
# You can changed the values above and rerun this program
# ---------------------------------------------------------------------------
import sys
import os
import copy
test_program = 'example/user/ill_condition.py'
if sys.argv[0] != 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')
# ------------------------------------------------------------------------
# Note that the a, t values are not used for this example
def example_db (file_name) :
def fun_chi(a, t) :
return ('prior_chi', 'prior_dchi', None)
# ----------------------------------------------------------------------
# age table:
age_list = [ 0.0, 100.0 ]
#
# time table:
time_list = [ 1990.0, 2010.0 ]
#
# integrand table:
integrand_table = [
{ 'name': 'mtexcess' }
]
#
# node table:
node_table = [ { 'name':'world', 'parent':'' } ]
#
# weight table:
weight_table = list()
#
# covariate table: empty
covariate_table = list()
#
# mulcov table: empty
mulcov_table = list()
#
# avgint table: empty
avgint_table = list()
#
# nslist_dict:
nslist_dict = dict()
# ----------------------------------------------------------------------
# data table:
data_table = list()
#
# values that are the same for all data rows
row = {
'meas_value': chi_bar ,
'density': 'gaussian',
'weight': '',
'hold_out': False,
'age_lower': age_list[0],
'age_upper': age_list[0],
'time_lower': 2000.,
'time_upper': 2000.,
'integrand': 'mtexcess',
'meas_std': sigma,
'node': 'world'
}
data_table.append(row)
# ----------------------------------------------------------------------
# prior_table
prior_table = [
{ # prior_chi
'name': 'prior_chi',
'density': 'uniform',
'lower': 0.0,
'upper': 1.0,
'mean': 0.0
},{ # prior_dchi
'name': 'prior_dchi',
'density': 'gaussian',
'lower': None,
'upper': None,
'mean': 0.0,
'std': delta
}
]
# ----------------------------------------------------------------------
# smooth table
name = 'smooth_chi'
fun = fun_chi
smooth_table = [
{'name':name, 'age_id':[0,1], 'time_id':[0], 'fun':fun }
]
# ----------------------------------------------------------------------
# rate table:
rate_table = [
{ 'name': 'chi',
'parent_smooth': 'smooth_chi',
'child_smooth': None
}
]
# ----------------------------------------------------------------------
# option_table
option_table = [
{ 'name':'rate_case', 'value':'iota_zero_rho_zero' },
{ 'name':'parent_node_name', 'value':'world' },
{ 'name':'random_seed', 'value':'0' },
{ 'name':'quasi_fixed', 'value':'false' },
{ 'name':'derivative_test_fixed', 'value':'second-order' },
{ 'name':'max_num_iter_fixed', 'value':'100' },
{ 'name':'print_level_fixed', 'value':'0' },
{ 'name':'tolerance_fixed', 'value':'1e-8' }
]
# ----------------------------------------------------------------------
# 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
)
# ----------------------------------------------------------------------
return
# ===========================================================================
# create example.db
file_name = 'example.db'
example_db(file_name)
program = '../../devel/dismod_at'
#
# init
dismod_at.system_command_prc([ program, file_name, 'init' ])
#
# fit fixed
dismod_at.system_command_prc([ program, file_name, 'fit', 'fixed' ])
#
# read database
connection = dismod_at.create_connection(
file_name, new = False, readonly = True
)
age_table = dismod_at.get_table_dict(connection, 'age')
rate_table = dismod_at.get_table_dict(connection, 'rate')
var_table = dismod_at.get_table_dict(connection, 'var')
fit_var_table = dismod_at.get_table_dict(connection, 'fit_var')
integrand_table = dismod_at.get_table_dict(connection, 'integrand')
#
for var_id in range( len(var_table) ) :
row = var_table[var_id]
rate_id = row['rate_id']
age_id = row['age_id']
time_id = row['time_id']
fit_value = fit_var_table[var_id]['fit_var_value']
assert row['var_type'] == 'rate'
assert rate_table[rate_id]['rate_name'] == 'chi'
assert time_id == 0
if age_id == 0 :
# this point is informed by the data so it fits
assert abs( 1.0 - fit_value / chi_bar ) < 1e-6
else :
# this point is only informed by prior which has a large
# standard deviation compared to the data
assert abs( 1.0 - fit_value / chi_bar ) > 1.0
#
# rescale
dismod_at.system_command_prc([ program, file_name, 'set', 'scale_var', 'fit_var' ])
#
# fit fixed
dismod_at.system_command_prc([ program, file_name, 'fit', 'fixed' ])
#
# read new fit values
fit_var_table = dismod_at.get_table_dict(connection, 'fit_var')
#
for var_id in range( len(var_table) ) :
age_id = row['age_id']
fit_value = fit_var_table[var_id]['fit_var_value']
if age_id == 0 :
# this point is informed by the data so it fits
assert abs( 1.0 - fit_value / chi_bar ) < 1e-10
else :
# this point is only informed by prior which has a large
# standard deviation compared to the data
assert abs( 1.0 - fit_value / chi_bar ) < 1e-5
# -----------------------------------------------------------------------------
print('ill_condition.py: OK')
# -----------------------------------------------------------------------------