<s> import unittest <EOL> import kabuki <EOL> import numpy as np <EOL> import pymc as pm <EOL> import math <EOL> import scipy as sc <EOL> import collections <EOL> from pprint import pprint <EOL> from numpy . random import randn <EOL> from numpy import array , sqrt <EOL> from nose import SkipTest <EOL> from pandas import DataFrame <EOL> from time import time <EOL> import unittest <EOL> from kabuki . utils import stochastic_from_dist <EOL> def multi_normal_like ( values , vec_mu , tau ) : <EOL> """<STR_LIT>""" <EOL> logp = <NUM_LIT:0> <EOL> for i in range ( len ( vec_mu ) ) : <EOL> logp += pm . normal_like ( values [ i , : ] , vec_mu [ i ] , tau ) <EOL> return logp <EOL> MN = stochastic_from_dist ( name = "<STR_LIT>" , logp = multi_normal_like ) <EOL> class TestStepMethods ( unittest . TestCase ) : <EOL> def __init__ ( self , * args , ** kwargs ) : <EOL> super ( TestStepMethods , self ) . __init__ ( * args , ** kwargs ) <EOL> self . uniform_lb = <NUM_LIT> <EOL> self . uniform_ub = <NUM_LIT> <EOL> def runTest ( self ) : <EOL> return <EOL> def assert_results ( self , node , true_value , true_mean , true_std = None , <EOL> mean_tol = <NUM_LIT:0.1> , std_tol = <NUM_LIT> ) : <EOL> """<STR_LIT>""" <EOL> pprint ( node . stats ( ) ) <EOL> lb = node . stats ( ) [ '<STR_LIT>' ] [ <NUM_LIT> ] <EOL> ub = node . stats ( ) [ '<STR_LIT>' ] [ <NUM_LIT> ] <EOL> if not ( lb < true_value < ub ) : <EOL> print ( "<STR_LIT>" ) <EOL> print ( "<STR_LIT>" , true_value ) <EOL> print ( "<STR_LIT>" , node . stats ( ) [ '<STR_LIT>' ] [ <NUM_LIT:50> ] ) <EOL> print ( "<STR_LIT>" , node . stats ( ) [ '<STR_LIT>' ] ) <EOL> print ( "<STR_LIT>" , true_mean ) <EOL> if true_std is not None : <EOL> print ( "<STR_LIT>" , true_std ) <EOL> print ( "<STR_LIT>" , node . stats ( ) [ '<STR_LIT>' ] ) <EOL> np . testing . assert_allclose ( node . stats ( ) [ '<STR_LIT>' ] , true_mean , rtol = mean_tol ) <EOL> if true_std is not None : <EOL> np . testing . assert_allclose ( node . stats ( ) [ '<STR_LIT>' ] , true_std , rtol = std_tol ) <EOL> def normal_normal ( self , add_shift , sigma_0 , sigma_beta , sigma_y , true_mu , <EOL> n_subjs , avg_samples , seed , use_metropolis ) : <EOL> """<STR_LIT>""" <EOL> np . random . seed ( seed ) <EOL> tau_0 = sigma_0 ** - <NUM_LIT:2> <EOL> mu_0 = <NUM_LIT:0.> <EOL> nodes , size , x_values = self . create_nodes_for_normal_normal ( add_shift , tau_0 , mu_0 , sigma_beta , <EOL> sigma_y , true_mu , n_subjs , avg_samples ) <EOL> mu = nodes [ '<STR_LIT>' ] <EOL> mm = pm . MCMC ( nodes ) <EOL> if use_metropolis : <EOL> mm . sample ( <NUM_LIT> , <NUM_LIT> ) <EOL> else : <EOL> mm . use_step_method ( kabuki . steps . kNormalNormal , mu ) <EOL> mm . sample ( <NUM_LIT> ) <EOL> total_n = sum ( size ) <EOL> tau = sigma_y ** - <NUM_LIT:2> <EOL> sum_obs = sum ( [ sum ( x . value . flatten ( ) ) for x in mm . observed_stochastics ] ) <EOL> if add_shift : <EOL> tmp = sum ( array ( size ) * x_values ) <EOL> else : <EOL> tmp = <NUM_LIT:0> <EOL> tau_prime = tau_0 + total_n * tau <EOL> mu_prime = ( tau * ( sum_obs - tmp ) + mu_0 * tau_0 ) / tau_prime <EOL> true_std = <NUM_LIT:1.> / np . sqrt ( tau_prime ) <EOL> self . assert_results ( mu , true_mu , mu_prime , true_std , mean_tol = <NUM_LIT:0.1> , std_tol = <NUM_LIT:0.1> ) <EOL> return mm , mu_prime , true_std <EOL> def create_nodes_for_normal_normal ( self , add_shift , tau_0 , mu_0 , sigma_beta , <EOL> sigma_y , true_mu , n_subjs , avg_samples ) : <EOL> """<STR_LIT>""" <EOL> mu = pm . Normal ( '<STR_LIT>' , mu_0 , tau_0 ) <EOL> nodes = { '<STR_LIT>' : mu } <EOL> size = [ None ] * n_subjs <EOL> x_values = [ None ] * n_subjs <EOL> if add_shift : <EOL> b = [ ] <EOL> else : <EOL> b = None <EOL> for i in range ( n_subjs ) : <EOL> size [ i ] = int ( max ( <NUM_LIT:1> , avg_samples + randn ( ) * <NUM_LIT:10> ) ) <EOL> if add_shift : <EOL> x_values [ i ] = randn ( ) * sigma_beta <EOL> value = randn ( size [ i ] ) * sigma_y + true_mu + x_values [ i ] <EOL> x = pm . Lambda ( '<STR_LIT>' % i , lambda x = x_values [ i ] : x ) <EOL> y = pm . Normal ( '<STR_LIT>' % i , mu + x , sigma_y ** - <NUM_LIT:2> , value = value , observed = True ) <EOL> nodes [ '<STR_LIT>' % i ] = x <EOL> b . append ( x ) <EOL> else : <EOL> value = randn ( size [ i ] ) * sigma_y + true_mu <EOL> y = pm . Normal ( '<STR_LIT>' % i , mu , sigma_y ** - <NUM_LIT:2> , value = value , observed = True ) <EOL> nodes [ '<STR_LIT>' % i ] = y <EOL> return nodes , size , x_values <EOL> def normal_normal_bundle ( self , use_metropolis ) : <EOL> """<STR_LIT>""" <EOL> self . normal_normal ( add_shift = True , sigma_0 = <NUM_LIT> , sigma_beta = <NUM_LIT> , sigma_y = <NUM_LIT> , <EOL> true_mu = - <NUM_LIT> , n_subjs = <NUM_LIT:1> , avg_samples = <NUM_LIT:100> , seed = <NUM_LIT:1> , use_metropolis = use_metropolis ) <EOL> self . normal_normal ( add_shift = True , sigma_0 = <NUM_LIT> , sigma_beta = <NUM_LIT> , sigma_y = <NUM_LIT:2> , <EOL> true_mu = - <NUM_LIT> , n_subjs = <NUM_LIT:2> , avg_samples = <NUM_LIT:10> , seed = <NUM_LIT:2> , use_metropolis = use_metropolis ) <EOL> self . normal_normal ( add_shift = True , sigma_0 = <NUM_LIT> , sigma_beta = <NUM_LIT:1.> , sigma_y = <NUM_LIT> , <EOL> true_mu = - <NUM_LIT:1.> , n_subjs = <NUM_LIT:3> , avg_samples = <NUM_LIT:10> , seed = <NUM_LIT:3> , use_metropolis = use_metropolis ) <EOL> self . normal_normal ( add_shift = False , sigma_0 = <NUM_LIT:1.> , sigma_beta = <NUM_LIT:0.5> , sigma_y = <NUM_LIT:0.5> , <EOL> true_mu = - <NUM_LIT> , n_subjs = <NUM_LIT:4> , avg_samples = <NUM_LIT:50> , seed = <NUM_LIT:4> , use_metropolis = use_metropolis ) <EOL> self . normal_normal ( add_shift = False , sigma_0 = <NUM_LIT> , sigma_beta = <NUM_LIT> , sigma_y = <NUM_LIT> , <EOL> true_mu = - <NUM_LIT> , n_subjs = <NUM_LIT:5> , avg_samples = <NUM_LIT:50> , seed = <NUM_LIT:5> , use_metropolis = use_metropolis ) <EOL> self . normal_normal ( add_shift = False , sigma_0 = <NUM_LIT> , sigma_beta = <NUM_LIT> , sigma_y = <NUM_LIT> , <EOL> true_mu = <NUM_LIT> , n_subjs = <NUM_LIT:6> , avg_samples = <NUM_LIT:30> , seed = <NUM_LIT:6> , use_metropolis = use_metropolis ) <EOL> def test_normal_normal_solution ( self ) : <EOL> """<STR_LIT>""" <EOL> self . normal_normal_bundle ( use_metropolis = True ) <EOL> def test_kNormalNormal ( self ) : <EOL> """<STR_LIT>""" <EOL> self . normal_normal_bundle ( use_metropolis = False ) <EOL> def create_nodes_for_PriorNormalstd ( self , n_subjs , sigma_0 , mu_0 , prior , ** kwargs ) : <EOL> """<STR_LIT>""" <EOL> if prior is pm . Uniform : <EOL> sigma = pm . Uniform ( '<STR_LIT>' , self . uniform_lb , self . uniform_ub , value = <NUM_LIT:1.> ) <EOL> elif prior is kabuki . utils . HalfCauchy : <EOL> sigma = kabuki . utils . HalfCauchy ( '<STR_LIT>' , ** kwargs ) <EOL> x_values = [ ] <EOL> nodes = { '<STR_LIT>' : sigma } <EOL> if not isinstance ( mu_0 , collections . Sequence ) : <EOL> mu_0 = [ mu_0 ] <EOL> n_conds = len ( mu_0 ) <EOL> for i_cond in range ( n_conds ) : <EOL> for i in range ( n_subjs ) : <EOL> x_value = randn ( ) * sigma_0 + mu_0 [ i_cond ] <EOL> x_values . append ( x_value ) <EOL> x = pm . Normal ( '<STR_LIT>' % i , mu_0 [ i_cond ] , sigma ** - <NUM_LIT:2> , value = x_value , observed = True ) <EOL> nodes [ '<STR_LIT>' % ( i , i_cond ) ] = x <EOL> return nodes , x_values <EOL> def uniform_normalstd ( self , sigma_0 , mu_0 , n_subjs , seed , use_metropolis ) : <EOL> """<STR_LIT>""" <EOL> np . random . seed ( seed ) <EOL> nodes , x_values = self . create_nodes_for_PriorNormalstd ( n_subjs , sigma_0 , mu_0 , prior = pm . Uniform ) <EOL> sigma = nodes [ '<STR_LIT>' ] <EOL> mm = pm . MCMC ( nodes ) <EOL> if use_metropolis : <EOL> mm . sample ( <NUM_LIT> , <NUM_LIT> ) <EOL> else : <EOL> mm . use_step_method ( kabuki . steps . UniformPriorNormalstd , sigma ) <EOL> mm . sample ( <NUM_LIT> ) <EOL> alpha = ( n_subjs - <NUM_LIT:1> ) / <NUM_LIT> <EOL> beta = sum ( [ ( x - mu_0 ) ** <NUM_LIT:2> for x in x_values ] ) / <NUM_LIT> <EOL> true_mean = math . gamma ( alpha - <NUM_LIT:0.5> ) / math . gamma ( alpha ) * np . sqrt ( beta ) <EOL> anal_var = beta / ( alpha - <NUM_LIT:1> ) - true_mean ** <NUM_LIT:2> <EOL> true_std = np . sqrt ( anal_var ) <EOL> self . assert_results ( sigma , sigma_0 , true_mean , true_std ) <EOL> return mm <EOL> def uniform_normalstd_multiple_conds_with_shared_sigma ( self , sigma_0 , mu_0 , n_subjs , seed , use_metropolis ) : <EOL> """<STR_LIT>""" <EOL> np . random . seed ( seed ) <EOL> n_conds = len ( mu_0 ) <EOL> nodes , x_values = self . create_nodes_for_PriorNormalstd ( n_subjs , sigma_0 , mu_0 , prior = pm . Uniform ) <EOL> sigma = nodes [ '<STR_LIT>' ] <EOL> mm = pm . MCMC ( nodes ) <EOL> if use_metropolis : <EOL> mm . sample ( <NUM_LIT> , <NUM_LIT> ) <EOL> else : <EOL> mm . use_step_method ( kabuki . steps . UniformPriorNormalstd , sigma ) <EOL> mm . sample ( <NUM_LIT> ) <EOL> alpha = ( n_subjs * n_conds - <NUM_LIT:1> ) / <NUM_LIT> <EOL> beta = <NUM_LIT:0> <EOL> for i_cond in range ( n_conds ) : <EOL> cur_x_values = x_values [ i_cond * n_subjs : ( i_cond + <NUM_LIT:1> ) * n_subjs ] <EOL> beta += sum ( [ ( x - mu_0 [ i_cond ] ) ** <NUM_LIT:2> for x in cur_x_values ] ) / <NUM_LIT> <EOL> true_mean = math . gamma ( alpha - <NUM_LIT:0.5> ) / math . gamma ( alpha ) * np . sqrt ( beta ) <EOL> anal_var = beta / ( alpha - <NUM_LIT:1> ) - true_mean ** <NUM_LIT:2> <EOL> true_std = np . sqrt ( anal_var ) <EOL> self . assert_results ( sigma , sigma_0 , true_mean , true_std ) <EOL> return mm <EOL> def test_uniform_normalstd_numerical_solution ( self ) : <EOL> """<STR_LIT>""" <EOL> self . uniform_normalstd ( sigma_0 = <NUM_LIT:0.5> , mu_0 = <NUM_LIT:0> , n_subjs = <NUM_LIT:8> , seed = <NUM_LIT:1> , use_metropolis = True ) <EOL> self . uniform_normalstd ( sigma_0 = <NUM_LIT> , mu_0 = - <NUM_LIT:100> , n_subjs = <NUM_LIT:4> , seed = <NUM_LIT:2> , use_metropolis = True ) <EOL> self . uniform_normalstd ( sigma_0 = <NUM_LIT> , mu_0 = <NUM_LIT:2> , n_subjs = <NUM_LIT:5> , seed = <NUM_LIT:3> , use_metropolis = True ) <EOL> self . uniform_normalstd ( sigma_0 = <NUM_LIT> , mu_0 = - <NUM_LIT:4> , n_subjs = <NUM_LIT:7> , seed = <NUM_LIT:4> , use_metropolis = True ) <EOL> def test_UniformNormalstd_step_method ( self ) : <EOL> """<STR_LIT>""" <EOL> self . uniform_normalstd ( sigma_0 = <NUM_LIT:0.5> , mu_0 = <NUM_LIT:0> , n_subjs = <NUM_LIT:8> , seed = <NUM_LIT:1> , use_metropolis = False ) <EOL> self . uniform_normalstd ( sigma_0 = <NUM_LIT> , mu_0 = - <NUM_LIT:100> , n_subjs = <NUM_LIT:4> , seed = <NUM_LIT:2> , use_metropolis = False ) <EOL> self . uniform_normalstd ( sigma_0 = <NUM_LIT> , mu_0 = <NUM_LIT:2> , n_subjs = <NUM_LIT:5> , seed = <NUM_LIT:3> , use_metropolis = False ) <EOL> self . uniform_normalstd ( sigma_0 = <NUM_LIT> , mu_0 = - <NUM_LIT:4> , n_subjs = <NUM_LIT:7> , seed = <NUM_LIT:4> , use_metropolis = False ) <EOL> self . uniform_normalstd ( sigma_0 = <NUM_LIT> , mu_0 = <NUM_LIT:10> , n_subjs = <NUM_LIT:4> , seed = <NUM_LIT:5> , use_metropolis = False ) <EOL> def test_uniform_normalstd_with_multiple_condition_numerical_solution ( self ) : <EOL> """<STR_LIT>""" <EOL> self . uniform_normalstd_multiple_conds_with_shared_sigma ( sigma_0 = <NUM_LIT:0.5> , mu_0 = ( <NUM_LIT:0> , <NUM_LIT:2> , <NUM_LIT:10> ) , n_subjs = <NUM_LIT:8> , seed = <NUM_LIT:1> , use_metropolis = True ) <EOL> self . uniform_normalstd_multiple_conds_with_shared_sigma ( sigma_0 = <NUM_LIT> , mu_0 = ( - <NUM_LIT:100> , <NUM_LIT:3> ) , n_subjs = <NUM_LIT:4> , seed = <NUM_LIT:2> , use_metropolis = True ) <EOL> self . uniform_normalstd_multiple_conds_with_shared_sigma ( sigma_0 = <NUM_LIT> , mu_0 = ( <NUM_LIT:1> , <NUM_LIT:2> ) , n_subjs = <NUM_LIT:5> , seed = <NUM_LIT:3> , use_metropolis = True ) <EOL> self . uniform_normalstd_multiple_conds_with_shared_sigma ( sigma_0 = <NUM_LIT:0.5> , mu_0 = ( - <NUM_LIT:4> , - <NUM_LIT:3> , <NUM_LIT:2> , <NUM_LIT:1> , <NUM_LIT:0> ) , n_subjs = <NUM_LIT:7> , seed = <NUM_LIT:4> , use_metropolis = True ) <EOL> def test_UniformNormalstd_step_method_with_multiple_condition ( self ) : <EOL> """<STR_LIT>""" <EOL> self . uniform_normalstd_multiple_conds_with_shared_sigma ( sigma_0 = <NUM_LIT:0.5> , mu_0 = ( <NUM_LIT:0> , <NUM_LIT:2> , <NUM_LIT:10> ) , n_subjs = <NUM_LIT:8> , seed = <NUM_LIT:1> , use_metropolis = False ) <EOL> self . uniform_normalstd_multiple_conds_with_shared_sigma ( sigma_0 = <NUM_LIT> , mu_0 = ( - <NUM_LIT:100> , <NUM_LIT:3> ) , n_subjs = <NUM_LIT:4> , seed = <NUM_LIT:2> , use_metropolis = False ) <EOL> self . uniform_normalstd_multiple_conds_with_shared_sigma ( sigma_0 = <NUM_LIT> , mu_0 = ( <NUM_LIT:1> , <NUM_LIT:2> ) , n_subjs = <NUM_LIT:5> , seed = <NUM_LIT:3> , use_metropolis = False ) <EOL> self . uniform_normalstd_multiple_conds_with_shared_sigma ( sigma_0 = <NUM_LIT:0.5> , mu_0 = [ - <NUM_LIT:4> ] , n_subjs = <NUM_LIT:7> , seed = <NUM_LIT:4> , use_metropolis = False ) <EOL> def numerical_solution ( self , defective_posterior , lb , ub ) : <EOL> """<STR_LIT>""" <EOL> norm_factor = sc . integrate . quad ( defective_posterior , lb , ub ) [ <NUM_LIT:0> ] <EOL> moment = lambda x , n = <NUM_LIT:1> : defective_posterior ( x ) * ( x ** n ) / norm_factor <EOL> true_mean = sc . integrate . quad ( moment , lb , ub , args = ( <NUM_LIT:1> ) ) [ <NUM_LIT:0> ] <EOL> m2 = sc . integrate . quad ( moment , lb , ub , args = ( <NUM_LIT:2> ) ) [ <NUM_LIT:0> ] <EOL> anal_var = m2 - true_mean ** <NUM_LIT:2> <EOL> true_std = np . sqrt ( anal_var ) <EOL> return true_mean , true_std <EOL> def half_cauchy_normal_std ( self , sigma_0 = <NUM_LIT:1.> , mu_0 = <NUM_LIT:0.> , S = <NUM_LIT:10> , n_subjs = <NUM_LIT:8> , seed = <NUM_LIT:1> , <EOL> use_metropolis = False ) : <EOL> """<STR_LIT>""" <EOL> np . random . seed ( seed ) <EOL> nodes , x_values = self . create_nodes_for_PriorNormalstd ( n_subjs , sigma_0 , mu_0 , <EOL> prior = kabuki . utils . HalfCauchy , <EOL> S = S , value = <NUM_LIT:1> ) <EOL> sigma = nodes [ '<STR_LIT>' ] <EOL> mm = pm . MCMC ( nodes ) <EOL> if use_metropolis : <EOL> mm . sample ( <NUM_LIT> , <NUM_LIT> ) <EOL> else : <EOL> mm . use_step_method ( kabuki . steps . HCauchyPriorNormalstd , sigma ) <EOL> mm . sample ( <NUM_LIT:10> ) <EOL> beta = sum ( ( array ( x_values ) - mu_0 ) ** <NUM_LIT:2> ) / <NUM_LIT:2> <EOL> def defective_posterior ( x , n = n_subjs , beta = beta , S = S ) : <EOL> gammapdf = ( x ** <NUM_LIT:2> ) ** ( - n / <NUM_LIT> ) * np . exp ( - beta / ( x ** <NUM_LIT:2> ) ) <EOL> cauchy = S / ( x ** <NUM_LIT:2> + S ** <NUM_LIT:2> ) <EOL> return gammapdf * cauchy <EOL> true_mean , true_std = self . numerical_solution ( defective_posterior , <NUM_LIT:0> , np . inf ) <EOL> self . assert_results ( sigma , sigma_0 , true_mean , true_std ) <EOL> return mm <EOL> def half_cauchy_bundle ( self , use_metropolis ) : <EOL> self . half_cauchy_normal_std ( sigma_0 = <NUM_LIT> , mu_0 = <NUM_LIT:0> , n_subjs = <NUM_LIT:8> , seed = <NUM_LIT:1> , S = <NUM_LIT:5> , <EOL> use_metropolis = use_metropolis ) <EOL> self . half_cauchy_normal_std ( sigma_0 = <NUM_LIT:0.5> , mu_0 = - <NUM_LIT:100> , n_subjs = <NUM_LIT:4> , seed = <NUM_LIT:2> , S = <NUM_LIT:20> , <EOL> use_metropolis = use_metropolis ) <EOL> self . half_cauchy_normal_std ( sigma_0 = <NUM_LIT> , mu_0 = <NUM_LIT:2> , n_subjs = <NUM_LIT:5> , seed = <NUM_LIT:3> , S = <NUM_LIT:3> , <EOL> use_metropolis = use_metropolis ) <EOL> self . half_cauchy_normal_std ( sigma_0 = <NUM_LIT> , mu_0 = - <NUM_LIT:4> , n_subjs = <NUM_LIT:7> , seed = <NUM_LIT:4> , S = <NUM_LIT:10> , <EOL> use_metropolis = use_metropolis ) <EOL> self . half_cauchy_normal_std ( sigma_0 = <NUM_LIT> , mu_0 = <NUM_LIT:10> , n_subjs = <NUM_LIT:4> , seed = <NUM_LIT:5> , S = <NUM_LIT:15> , <EOL> use_metropolis = use_metropolis ) <EOL> def test_half_cauchy_numerical_solution ( self ) : <EOL> """<STR_LIT>""" <EOL> self . half_cauchy_bundle ( use_metropolis = True ) <EOL> def test_HCauchyNormalstd_step_method ( self ) : <EOL> """<STR_LIT>""" <EOL> raise SkipTest ( "<STR_LIT>" ) <EOL> self . half_cauchy_bundle ( use_metropolis = False ) <EOL> @ unittest . skip ( "<STR_LIT>" ) <EOL> def run_SPXcentered ( self , sigma_x , n_subjs , size , mu_value , mu_step_method , seed ) : <EOL> """<STR_LIT>""" <EOL> if np . isscalar ( mu_value ) : <EOL> n_conds = <NUM_LIT:1> <EOL> else : <EOL> n_conds = len ( mu_value ) <EOL> max_tries = <NUM_LIT:5> <EOL> iter = <NUM_LIT> <EOL> burnin = <NUM_LIT> <EOL> nodes , t_values = self . create_hierarchical_model ( sigma_x = sigma_x , n_subjs = n_subjs , size = size , <EOL> mu_value = mu_value , seed = seed ) <EOL> mcmc = pm . MCMC ( nodes ) <EOL> [ mcmc . use_step_method ( mu_step_method , node ) for node in nodes [ '<STR_LIT>' ] ] <EOL> nodes_spx , t_values = self . create_hierarchical_model ( sigma_x = sigma_x , n_subjs = n_subjs , size = size , <EOL> mu_value = mu_value , seed = seed ) <EOL> mcmc_spx = pm . MCMC ( nodes_spx ) <EOL> mcmc_spx . use_step_method ( kabuki . steps . SPXcentered , loc = nodes_spx [ '<STR_LIT>' ] , <EOL> scale = nodes_spx [ '<STR_LIT>' ] , <EOL> loc_step_method = mu_step_method ) <EOL> nodes_vpx , t_values = self . create_hierarchical_model ( sigma_x = sigma_x , n_subjs = n_subjs , size = size , <EOL> mu_value = mu_value , seed = seed , vec = True ) <EOL> mcmc_vpx = pm . MCMC ( nodes_vpx ) <EOL> mcmc_vpx . use_step_method ( kabuki . steps . SPXcentered , loc = nodes_vpx [ '<STR_LIT>' ] , <EOL> scale = nodes_vpx [ '<STR_LIT>' ] , <EOL> loc_step_method = mu_step_method ) <EOL> i_try = <NUM_LIT:0> <EOL> while i_try < max_tries : <EOL> print ( "<STR_LIT>" % ( i_try + <NUM_LIT:1> ) ) <EOL> i_t = time ( ) <EOL> mcmc_spx . sample ( iter , burnin ) <EOL> print ( "<STR_LIT>" % ( time ( ) - i_t ) ) <EOL> stats = dict ( [ ( '<STR_LIT>' % x , mcmc_spx . mu [ x ] . stats ( ) ) for x in range ( n_conds ) ] ) <EOL> i_t = time ( ) <EOL> mcmc_vpx . sample ( iter , burnin ) <EOL> print ( "<STR_LIT>" % ( time ( ) - i_t ) ) <EOL> stats . update ( dict ( [ ( '<STR_LIT>' % x , mcmc_vpx . mu [ x ] . stats ( ) ) for x in range ( n_conds ) ] ) ) <EOL> i_t = time ( ) <EOL> mcmc . sample ( iter , burnin ) <EOL> print ( "<STR_LIT>" % ( time ( ) - i_t ) ) <EOL> stats . update ( dict ( [ ( '<STR_LIT>' % x , mcmc . mu [ x ] . stats ( ) ) for x in range ( n_conds ) ] ) ) <EOL> df = DataFrame ( stats , index = [ '<STR_LIT>' , '<STR_LIT>' ] ) . T <EOL> df = df . rename ( columns = { '<STR_LIT>' : '<STR_LIT>' , '<STR_LIT>' : '<STR_LIT>' } ) <EOL> print ( df ) <EOL> try : <EOL> for i in range ( len ( df ) / <NUM_LIT:3> ) : <EOL> np . testing . assert_allclose ( df [ ( <NUM_LIT:3> * i + <NUM_LIT:0> ) : ( <NUM_LIT:3> * i + <NUM_LIT:1> ) ] , df [ ( <NUM_LIT:3> * i + <NUM_LIT:1> ) : ( <NUM_LIT:3> * i + <NUM_LIT:2> ) ] , atol = <NUM_LIT:0.1> , rtol = <NUM_LIT> ) <EOL> np . testing . assert_allclose ( df [ ( <NUM_LIT:3> * i + <NUM_LIT:1> ) : ( <NUM_LIT:3> * i + <NUM_LIT:2> ) ] , df [ ( <NUM_LIT:3> * i + <NUM_LIT:2> ) : ( <NUM_LIT:3> * i + <NUM_LIT:3> ) ] , atol = <NUM_LIT:0.1> , rtol = <NUM_LIT> ) <EOL> np . testing . assert_allclose ( df [ ( <NUM_LIT:3> * i + <NUM_LIT:2> ) : ( <NUM_LIT:3> * i + <NUM_LIT:3> ) ] , df [ ( <NUM_LIT:3> * i + <NUM_LIT:0> ) : ( <NUM_LIT:3> * i + <NUM_LIT:1> ) ] , atol = <NUM_LIT:0.1> , rtol = <NUM_LIT> ) <EOL> break <EOL> except AssertionError : <EOL> print ( "<STR_LIT>" ) <EOL> i_try += <NUM_LIT:1> <EOL> assert ( i_try < max_tries ) , "<STR_LIT>" <EOL> @ unittest . skip ( "<STR_LIT>" ) <EOL> def test_SPX ( self ) : <EOL> """<STR_LIT>""" <EOL> print ( "<STR_LIT>" ) <EOL> self . run_SPXcentered ( sigma_x = <NUM_LIT:1> , n_subjs = <NUM_LIT:5> , size = <NUM_LIT:100> , mu_value = <NUM_LIT:4> , <EOL> mu_step_method = kabuki . steps . kNormalNormal , seed = <NUM_LIT:1> ) <EOL> print ( "<STR_LIT>" ) <EOL> self . run_SPXcentered ( sigma_x = <NUM_LIT:1> , n_subjs = <NUM_LIT:5> , size = <NUM_LIT:10> , mu_value = ( <NUM_LIT:4> , <NUM_LIT:3> , <NUM_LIT:2> , <NUM_LIT:1> , <NUM_LIT:0> , <NUM_LIT:4> , <NUM_LIT:3> , <NUM_LIT:2> , <NUM_LIT:1> , <NUM_LIT:0> ) , <EOL> mu_step_method = kabuki . steps . kNormalNormal , seed = <NUM_LIT:1> ) <EOL> print ( "<STR_LIT>" ) <EOL> self . run_SPXcentered ( sigma_x = <NUM_LIT:0.5> , n_subjs = <NUM_LIT:5> , size = <NUM_LIT:10> , mu_value = ( <NUM_LIT:4> , <NUM_LIT:3> ) , <EOL> mu_step_method = kabuki . steps . kNormalNormal , seed = <NUM_LIT:1> ) <EOL> print ( "<STR_LIT>" ) <EOL> self . run_SPXcentered ( sigma_x = <NUM_LIT:0.1> , n_subjs = <NUM_LIT:5> , size = <NUM_LIT:10> , mu_value = ( <NUM_LIT:4> , <NUM_LIT:3> ) , <EOL> mu_step_method = kabuki . steps . kNormalNormal , seed = <NUM_LIT:1> ) <EOL> print ( "<STR_LIT>" ) <EOL> self . run_SPXcentered ( sigma_x = <NUM_LIT:1> , n_subjs = <NUM_LIT:5> , size = <NUM_LIT:10> , mu_value = list ( range ( <NUM_LIT:20> ) ) , <EOL> mu_step_method = kabuki . steps . kNormalNormal , seed = <NUM_LIT:1> ) <EOL> print ( "<STR_LIT>" ) <EOL> self . run_SPXcentered ( sigma_x = <NUM_LIT:0.1> , n_subjs = <NUM_LIT:5> , size = <NUM_LIT:10> , mu_value = list ( range ( <NUM_LIT:20> ) ) , <EOL> mu_step_method = kabuki . steps . kNormalNormal , seed = <NUM_LIT:1> ) <EOL> def create_hierarchical_model ( self , sigma_x = <NUM_LIT:1> , n_subjs = <NUM_LIT:5> , size = <NUM_LIT:100> , mu_value = <NUM_LIT:4> , seed = <NUM_LIT:1> , vec = False ) : <EOL> """<STR_LIT>""" <EOL> np . random . seed ( seed ) <EOL> if np . isscalar ( mu_value ) : <EOL> mu_value = [ mu_value ] <EOL> n_conds = len ( mu_value ) <EOL> mu = [ None ] * n_conds <EOL> subj_nodes = [ None ] * n_conds <EOL> data_nodes = [ None ] * n_conds <EOL> true_values = { } <EOL> true_values [ '<STR_LIT>' ] = sigma_x <EOL> true_values [ '<STR_LIT>' ] = mu_value <EOL> sigma = pm . Uniform ( '<STR_LIT>' , <NUM_LIT> , <NUM_LIT> , value = <NUM_LIT:1> ) <EOL> tau = sigma ** - <NUM_LIT:2> <EOL> for i_cond in range ( n_conds ) : <EOL> true_x = randn ( n_subjs ) * sigma_x + mu_value [ i_cond ] <EOL> value = np . random . randn ( n_subjs , size ) . T + true_x <EOL> value = value . T <EOL> print ( true_x ) <EOL> mu [ i_cond ] = pm . Normal ( '<STR_LIT>' % i_cond , <NUM_LIT:0> , <NUM_LIT> ** - <NUM_LIT:2> , value = <NUM_LIT:0> ) <EOL> if vec : <EOL> subj_nodes [ i_cond ] = pm . Normal ( '<STR_LIT>' % ( i_cond ) , mu [ i_cond ] , tau , size = n_subjs ) <EOL> data_nodes [ i_cond ] = MN ( '<STR_LIT>' % ( i_cond ) , vec_mu = subj_nodes [ i_cond ] , tau = <NUM_LIT:1> , value = value , observed = True ) <EOL> else : <EOL> subj_nodes [ i_cond ] = [ None ] * n_subjs <EOL> data_nodes [ i_cond ] = [ None ] * n_subjs <EOL> for i_subj in range ( n_subjs ) : <EOL> subj_nodes [ i_cond ] [ i_subj ] = pm . Normal ( '<STR_LIT>' % ( i_cond , i_subj ) , mu [ i_cond ] , tau ) <EOL> data_nodes [ i_cond ] [ i_subj ] = pm . Normal ( '<STR_LIT>' % ( i_cond , i_subj ) , <EOL> mu = subj_nodes [ i_cond ] [ i_subj ] , <EOL> tau = <NUM_LIT:1> , value = value [ i_subj , : ] , observed = True ) <EOL> nodes = { } <EOL> nodes [ '<STR_LIT:x>' ] = subj_nodes <EOL> nodes [ '<STR_LIT:y>' ] = data_nodes <EOL> nodes [ '<STR_LIT>' ] = mu <EOL> nodes [ '<STR_LIT>' ] = sigma <EOL> return nodes , true_values <EOL> def run_SliceStep ( self , sigma_x , n_subjs , size , mu_value , seed , left = None , max_tries = <NUM_LIT:5> ) : <EOL> if np . isscalar ( mu_value ) : <EOL> n_conds = <NUM_LIT:1> <EOL> else : <EOL> n_conds = len ( mu_value ) <EOL> iter = <NUM_LIT> <EOL> burnin = <NUM_LIT> <EOL> nodes , t_values = self . create_hierarchical_model ( sigma_x = sigma_x , n_subjs = n_subjs , size = size , <EOL> mu_value = mu_value , seed = seed ) <EOL> mcmc = pm . MCMC ( nodes ) <EOL> [ mcmc . use_step_method ( kabuki . steps . kNormalNormal , node ) for node in nodes [ '<STR_LIT>' ] ] <EOL> nodes_s , t_values = self . create_hierarchical_model ( sigma_x = sigma_x , n_subjs = n_subjs , size = size , <EOL> mu_value = mu_value , seed = seed ) <EOL> mcmc_s = pm . MCMC ( nodes_s ) <EOL> [ mcmc_s . use_step_method ( kabuki . steps . kNormalNormal , node ) for node in nodes_s [ '<STR_LIT>' ] ] <EOL> mcmc_s . use_step_method ( kabuki . steps . SliceStep , nodes_s [ '<STR_LIT>' ] , width = <NUM_LIT:3> , left = left ) <EOL> i_try = <NUM_LIT:0> <EOL> stats = { } <EOL> while i_try < max_tries : <EOL> print ( "<STR_LIT>" % ( i_try + <NUM_LIT:1> ) ) <EOL> i_t = time ( ) <EOL> mcmc_s . sample ( iter , burnin ) <EOL> print ( "<STR_LIT>" % ( time ( ) - i_t ) ) <EOL> stats . update ( dict ( [ ( '<STR_LIT>' % x , mcmc_s . mu [ x ] . stats ( ) ) for x in range ( n_conds ) ] ) ) <EOL> i_t = time ( ) <EOL> mcmc . sample ( iter , burnin ) <EOL> print ( "<STR_LIT>" % ( time ( ) - i_t ) ) <EOL> stats . update ( dict ( [ ( '<STR_LIT>' % x , mcmc . mu [ x ] . stats ( ) ) for x in range ( n_conds ) ] ) ) <EOL> df = DataFrame ( stats , index = [ '<STR_LIT>' , '<STR_LIT>' ] ) . T <EOL> df = df . rename ( columns = { '<STR_LIT>' : '<STR_LIT>' , '<STR_LIT>' : '<STR_LIT>' } ) <EOL> print ( df ) <EOL> try : <EOL> for i in range ( len ( df ) / <NUM_LIT:2> ) : <EOL> np . testing . assert_allclose ( df [ ( <NUM_LIT:2> * i + <NUM_LIT:0> ) : ( <NUM_LIT:2> * i + <NUM_LIT:1> ) ] , df [ ( <NUM_LIT:2> * i + <NUM_LIT:1> ) : ( <NUM_LIT:2> * i + <NUM_LIT:2> ) ] , atol = <NUM_LIT:0.1> , rtol = <NUM_LIT> ) <EOL> break <EOL> except AssertionError : <EOL> print ( "<STR_LIT>" ) <EOL> print ( df [ ( <NUM_LIT:2> * i ) : ( <NUM_LIT:2> * ( i + <NUM_LIT:1> ) ) ] ) <EOL> print ( "<STR_LIT>" ) <EOL> i_try += <NUM_LIT:1> <EOL> assert ( i_try < max_tries ) , "<STR_LIT>" <EOL> return mcmc , mcmc_s <EOL> def test_SliceStep ( self ) : <EOL> """<STR_LIT>""" <EOL> print ( "<STR_LIT>" ) <EOL> self . run_SliceStep ( sigma_x = <NUM_LIT:1> , n_subjs = <NUM_LIT:5> , size = <NUM_LIT:100> , mu_value = <NUM_LIT:4> , seed = <NUM_LIT:1> ) <EOL> print ( "<STR_LIT>" ) <EOL> self . run_SliceStep ( sigma_x = <NUM_LIT:1> , n_subjs = <NUM_LIT:5> , size = <NUM_LIT:10> , mu_value = list ( range ( <NUM_LIT:10> ) ) , seed = <NUM_LIT:1> ) </s>
