1from ..generic import Base_Generic 2from ..generic.sample import logp_rho_prec, logp_lambda_prec 3from ...utils import se_covariance, se_precision 4from ... import verify 5import numpy as np 6 7class Base_SESE(Base_Generic): 8 """ 9 All arguments to this function ar documented in SESE 10 """ 11 def __init__(self, Y, X, W, M, Delta, 12 n_samples=1000, n_jobs=1, 13 extra_traced_params = None, 14 priors=None, 15 configs=None, 16 starting_values=None, 17 truncation=None): 18 super(Base_SESE, self).__init__(Y, X, W, M, Delta, 19 n_samples=0, n_jobs=n_jobs, 20 extra_traced_params=extra_traced_params, 21 priors=priors, 22 configs=configs, 23 starting_values=starting_values, 24 truncation=truncation) 25 self.state.Psi_1 = se_covariance 26 self.state.Psi_2 = se_covariance 27 self.state.Psi_1i = se_precision 28 self.state.Psi_2i = se_precision 29 30 self.configs.Rho.logp = logp_rho_prec 31 self.configs.Lambda.logp = logp_lambda_prec 32 33 try: 34 self.sample(n_samples, n_jobs =n_jobs) 35 except (np.linalg.LinAlgError, ValueError) as e: 36 warn('Encountered the following LinAlgError. ' 37 'Model will return for debugging purposes. \n {}'.format(e)) 38 39class SESE(Base_SESE): 40 """ 41 This is a class that provides the generic structures required for the two-level variance components model with simultaneous autregressive errors in both levels. 42 43 Formally, this is the following distributional model for Y: 44 45 Y ~ N(X * Beta, Phi_1(Rho, Sigma2) + Delta Phi_2(Lambda, Tau2) Delta') 46 47 Where Delta is the dummy variable matrix, Rho, Sigma2 are response-level autoregressive and scale components of a spatial autoregressive covariance function, Phi_1. Lambda and Tau2 are regional-level components for spatial autoregressive covariance function, Phi_2. In this case, Phi_1 and Phi_2 are Simultaneously-autoregressive errors over weights matrices W,M: 48 49 Phi_1(Rho, Sigma2) = [(I - Rho W)'(I - Rho W)]^{-1} * Sigma2 50 Phi_2(Lambda, Tau2) = [(I - Lambda M)'(I - Lambda M)]^{-1} * Tau2 51 52 Arguments 53 ---------- 54 Y : numpy.ndarray 55 The (n,1) array containing the response to be modeled 56 X : numpy.ndarray 57 The (n,p) array containing covariates used to predict the response, Y. 58 W : pysal.weights.W 59 a spatial weights object for the n observations 60 M : pysal.weights.W 61 a spatial weights object for the J regions. 62 Z : numpy.ndarray 63 The (J,p') array of region-level covariates used to predict the response, Y. 64 Delta : numpy.ndarray 65 The (n,J) dummy variable matrix, relating observation i to region j. Required if membership is not passed. 66 membership: numpy.ndarray 67 The (n,) vector of labels assigning each observation to a region. Required if Delta is not passed. 68 transform : string 69 weights transform to use in modeling. 70 verbose : bool 71 whether to provide verbose output about sampling progress 72 n_samples : int 73 the number of samples to draw. If zero, the model will only be initialized, and later sampling can be done using model.sample 74 n_jobs : int 75 the number of parallel chains to run. Defaults to 1. 76 extra_traced_param : list of strings 77 extra parameters in the trace to record. 78 center : bool 79 Whether to center the input data so that its mean is zero. Occasionally improves the performance of the sampler. 80 scale : bool 81 Whether to rescale the input data so that its variance is one. May dramatically improve the performance of the sampler, at the cost of interpretability of the coefficients. 82 priors : dictionary 83 A dictionary used to configure the priors for the model. This may include the following keys: 84 Betas_cov0 : prior covariance for Beta parameter vector 85 (default: I*100) 86 Betas_mean0 : prior mean for Beta parameters 87 (default: 0) 88 Sigma2_a0 : prior shape parameter for inverse gamma prior on response-level variance 89 (default: .001) 90 Sigma2_b0 : prior scale parameter for inverse gamma prior on response-level variance 91 (default: .001) 92 Tau2_a0 : prior shape parameter for inverse gamma prior on regional-level variance 93 (default: .001) 94 Tau2_b0 : prior scale parameter for inverse gamma prior on regional-level variance 95 (default: .001) 96 Log_Lambda0 : prior on Lambda, the region-level autoregressive parameter. Must be a callable function that takes a single argument and returns a single value providing the log prior likelihood. 97 (default: priors.constant) 98 Log_Rho0 : prior on Rho, the response-level autoregressive paraameter. Must be a callable function that takes a single argument and returns a single value providing the log prior likelihood. 99 (default: priors.constant) 100 starting_value : dictionary 101 A dictionary containing the starting values of the sampler. If n_jobs > 1, these starting values are perturbed randomly to induce overdispersion. 102 103 This dicutionary may include the following keys: 104 Betas : starting vector of Beta parameters. 105 (default: np.zeros(p,1)) 106 Alphas : starting vector of Alpha variance components. 107 (default: np.zeros(J,1)) 108 Sigma2 : starting value of response-level variance 109 (default: 4) 110 Tau2 : starting value of region-level varaince 111 (default: 4) 112 Rho : starting value of response-level spatial autoregressive parameter 113 (default: -1/n) 114 Lambda : starting value of region-level spatial autoregressive parameter 115 (default -1/J) 116 config : dictionary 117 A dictionary containing the configuration values for the non-Gibbs samplers for the spatial parameters. May be specified for each parameter individually (if both are in the model), or may be specified implicitly for the relevant parameter to the model. Keys may include: 118 Rho_method : string specifying whether "met" or "slice" sampling should be used for the response-level effects 119 Rho_configs : configuration options for the sampler for the response-level effects that will be used in the step method. 120 Lambda_method : string specifying whether 'met' or 'slice' sampling should be used for the region-level effects 121 Lambda_configs : configuration options for the sampler for the region-level effects that will be used in the step method. 122 123 For options that can be in Lambda/Rho_configs, see: 124 spvcm.steps.Slice, spvcm.steps.Metropolis 125 truncation : dictionary 126 A dictionary containing the configuration values for the maximum and minimum allowable Lambda and Rho parameters. If these are not provided, the support is row-standardized by default, and the minimal eigenvalue computed for the lower bound on the parameters. *only* the single minimum eigenvalue is computed, so this is still rather efficient for large matrices. Keys may include: 127 Rho_min : minimum value allowed for response-level autoregressive coefficient 128 Rho_max : maximum value allowed for response-level autoregressive coefficient 129 Lambda_min : minimum value allowed for region-level autoregressive coefficient 130 Lambda_max : maximum value allowed for region-level autoregressive coefficient. 131 """ 132 def __init__(self, Y, X, W, M, Z=None, Delta=None, membership=None, 133 #data options 134 transform ='r', verbose=False, 135 n_samples=1000, n_jobs=1, 136 extra_traced_params = None, 137 priors=None, 138 configs=None, 139 starting_values=None, 140 truncation=None, 141 center=False, 142 scale=False): 143 if X is None: 144 X = np.ones_like(Y) 145 center=False 146 scale=False 147 W,M = verify.weights(W, M, transform=transform) 148 self.M = M 149 N,_ = X.shape 150 J = M.n 151 Mmat = M.sparse 152 Wmat = W.sparse 153 154 Delta, membership = verify.Delta_members(Delta, membership, N, J) 155 if Z is not None: 156 Z = Delta.dot(Z) 157 X = np.hstack((X,Z)) 158 if center: 159 X = verify.center(X) 160 if scale: 161 X = verify.scale(X) 162 163 164 X = verify.covariates(X) 165 166 self._verbose = verbose 167 168 super(SESE, self).__init__(Y, X, Wmat, Mmat, Delta, 169 n_samples=n_samples, 170 n_jobs = n_jobs, 171 extra_traced_params=extra_traced_params, 172 priors=priors, 173 configs=configs, 174 starting_values=starting_values, 175 truncation=truncation) 176