1from __future__ import division 2 3import copy 4import scipy.stats as stats 5import numpy as np 6import numpy.linalg as la 7 8from ...both_levels.generic import Base_Generic 9from ...both_levels.generic.model import SAMPLERS as generic_parameters 10from ... import verify 11from ...utils import se_covariance, se_precision, ind_covariance, splogdet, chol_mvn 12from spreg.utils import spdot 13 14 15 16SAMPLERS = ['Alphas', 'Betas', 'Sigma2', 'Tau2', 'Lambda'] 17 18class Base_Upper_SE(Base_Generic): 19 """ 20 All arguments are documented in Upper_SE 21 """ 22 def __init__(self, Y, X, M, Delta, n_samples=1000, n_jobs=1, 23 extra_traced_params = None, 24 priors=None, 25 starting_values=None, 26 configs=None, 27 truncation=None): 28 W = np.eye((Delta.shape[0])) 29 super(Base_Upper_SE, self).__init__(Y, X, W, M, Delta, 30 n_samples=0, n_jobs=n_jobs, 31 extra_traced_params=None, 32 priors=priors, 33 configs=configs, 34 starting_values=starting_values, 35 truncation=truncation) 36 self.state.Psi_1 = ind_covariance 37 self.state.Psi_1i = ind_covariance 38 self.state.Psi_2 = se_covariance 39 self.state.Psi_2i = se_precision 40 original_traced = copy.deepcopy(self.traced_params) 41 to_drop = [k for k in original_traced if (k not in SAMPLERS and k in generic_parameters)] 42 self.traced_params = copy.deepcopy(SAMPLERS) 43 for param in to_drop: 44 for i, _ in enumerate(self.trace.chains): 45 del self.trace.chains[i][param] 46 extra_traced_params = [] if extra_traced_params is None else extra_traced_params 47 for extra in extra_traced_params: 48 self.traced_params.append(extra) 49 for i, chain in enumerate(self.trace.chains): 50 chain[extra] = [] 51 52 if n_samples > 0: 53 try: 54 self.sample(n_samples, n_jobs=n_jobs) 55 except (np.linalg.LinAlgError, ValueError) as e: 56 Warn('Encountered the following LinAlgError. ' 57 'Model will return for debugging. \n {}'.format(e)) 58 59 def _iteration(self): 60 st = self.state 61 62 ### Sample the Beta conditional posterior 63 ### P(beta | . ) \propto L(Y|.) \dot P(\beta) 64 ### is 65 ### N(Sb, S) where 66 ### S = (X' Sigma^{-1}_Y X + S_0^{-1})^{-1} 67 ### b = X' Sigma^{-1}_Y (Y - Delta Alphas) + S^{-1}\mu_0 68 covm_update = st.X.T.dot(st.X) / st.Sigma2 69 covm_update += st.Betas_cov0i 70 covm_update = la.inv(covm_update) 71 72 resids = st.Y - st.Delta.dot(st.Alphas) 73 XtSresids = st.X.T.dot(resids) / st.Sigma2 74 mean_update = XtSresids + st.Betas_cov0i.dot(st.Betas_mean0) 75 mean_update = np.dot(covm_update, mean_update) 76 st.Betas = chol_mvn(mean_update, covm_update) 77 st.XBetas = np.dot(st.X, st.Betas) 78 79 ### Sample the Random Effect conditional posterior 80 ### P( Alpha | . ) \propto L(Y|.) \dot P(Alpha | \lambda, Tau2) 81 ### \dot P(Tau2) \dot P(\lambda) 82 ### is 83 ### N(Sb, S) 84 ### Where 85 ### S = (Delta'Sigma_Y^{-1}Delta + Sigma_Alpha^{-1})^{-1} 86 ### b = (Delta'Sigma_Y^{-1}(Y - X\beta) + 0) 87 covm_update = st.Delta.T.dot(st.Delta) / st.Sigma2 88 covm_update += st.PsiLambdai / st.Tau2 89 covm_update = np.asarray(covm_update) 90 covm_update = la.inv(covm_update) 91 92 resids = st.Y - st.XBetas 93 mean_update = st.Delta.T.dot(resids) / st.Sigma2 94 mean_update = np.dot(covm_update, mean_update) 95 mean_update = np.asarray(mean_update) 96 st.Alphas = chol_mvn(mean_update, covm_update) 97 st.DeltaAlphas = np.dot(st.Delta, st.Alphas) 98 99 ### Sample the Random Effect aspatial variance parameter 100 ### P(Tau2 | .) \propto L(Y|.) \dot P(\Alpha | \lambda, Tau2) 101 ### \dot P(Tau2) \dot P(\lambda) 102 ### is 103 ### IG(J/2 + a0, u'(\Psi(\lambda))^{-1}u * .5 + b0) 104 bn = spdot(st.Alphas.T, spdot(st.PsiLambdai, st.Alphas)) * .5 + st.Tau2_b0 105 st.Tau2 = stats.invgamma.rvs(st.Tau2_an, scale=bn) 106 107 ### Sample the response aspatial variance parameter 108 ### P(Sigma2 | . ) \propto L(Y | .) \dot P(Sigma2) 109 ### is 110 ### IG(N/2 + a0, eta'Psi(\rho)^{-1}eta * .5 + b0) 111 ### Where eta is the linear predictor, Y - X\beta + \DeltaAlphas 112 eta = st.Y - st.XBetas - st.DeltaAlphas 113 bn = eta.T.dot(eta) * .5 + st.Sigma2_b0 114 st.Sigma2 = stats.invgamma.rvs(st.Sigma2_an, scale=bn) 115 116 ### P(Psi(\rho) | . ) \propto L(Y | .) \dot P(\rho) 117 ### is 118 ### |Psi(rho)|^{-1/2} exp(1/2(eta'Psi(rho)^{-1}eta * Sigma2^{-1})) * 1/(emax-emin) 119 st.Lambda = self.configs.Lambda(st) 120 st.PsiLambdai = st.Psi_2i(st.Lambda, st.M) 121 122class Upper_SE(Base_Upper_SE): 123 """ 124 This is a class that provides the generic structures required for the two-level variance components model with simultaneous autregressive errors in the region level and independent errors in the response level. 125 126 Formally, this is the following distributional model for Y: 127 128 Y ~ N(X * Beta, Phi_1(Rho, Sigma2) + Delta Phi_2(Lambda, Tau2) Delta') 129 130 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: 131 132 Phi_1(Rho, Sigma2) = I * Sigma2 133 Phi_2(Lambda, Tau2) = [(I - Lambda M)'(I - Lambda M)]^{-1} * Tau2 134 135 Arguments 136 ---------- 137 Y : numpy.ndarray 138 The (n,1) array containing the response to be modeled 139 X : numpy.ndarray 140 The (n,p) array containing covariates used to predict the response, Y. 141 W : pysal.weights.W 142 a spatial weights object for the n observations 143 M : pysal.weights.W 144 a spatial weights object for the J regions. 145 Z : numpy.ndarray 146 The (J,p') array of region-level covariates used to predict the response, Y. 147 Delta : numpy.ndarray 148 The (n,J) dummy variable matrix, relating observation i to region j. Required if membership is not passed. 149 membership: numpy.ndarray 150 The (n,) vector of labels assigning each observation to a region. Required if Delta is not passed. 151 transform : string 152 weights transform to use in modeling. 153 verbose : bool 154 whether to provide verbose output about sampling progress 155 n_samples : int 156 the number of samples to draw. If zero, the model will only be initialized, and later sampling can be done using model.sample 157 n_jobs : int 158 the number of parallel chains to run. Defaults to 1. 159 extra_traced_param : list of strings 160 extra parameters in the trace to record. 161 center : bool 162 Whether to center the input data so that its mean is zero. Occasionally improves the performance of the sampler. 163 scale : bool 164 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. 165 priors : dictionary 166 A dictionary used to configure the priors for the model. This may include the following keys: 167 Betas_cov0 : prior covariance for Beta parameter vector 168 (default: I*100) 169 Betas_mean0 : prior mean for Beta parameters 170 (default: 0) 171 Sigma2_a0 : prior shape parameter for inverse gamma prior on response-level variance 172 (default: .001) 173 Sigma2_b0 : prior scale parameter for inverse gamma prior on response-level variance 174 (default: .001) 175 Tau2_a0 : prior shape parameter for inverse gamma prior on regional-level variance 176 (default: .001) 177 Tau2_b0 : prior scale parameter for inverse gamma prior on regional-level variance 178 (default: .001) 179 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. 180 (default: priors.constant) 181 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. 182 (default: priors.constant) 183 starting_value : dictionary 184 A dictionary containing the starting values of the sampler. If n_jobs > 1, these starting values are perturbed randomly to induce overdispersion. 185 186 This dicutionary may include the following keys: 187 Betas : starting vector of Beta parameters. 188 (default: np.zeros(p,1)) 189 Alphas : starting vector of Alpha variance components. 190 (default: np.zeros(J,1)) 191 Sigma2 : starting value of response-level variance 192 (default: 4) 193 Tau2 : starting value of region-level varaince 194 (default: 4) 195 Rho : starting value of response-level spatial autoregressive parameter 196 (default: -1/n) 197 Lambda : starting value of region-level spatial autoregressive parameter 198 (default -1/J) 199 config : dictionary 200 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: 201 Rho_method : string specifying whether "met" or "slice" sampling should be used for the response-level effects 202 Rho_configs : configuration options for the sampler for the response-level effects that will be used in the step method. 203 Lambda_method : string specifying whether 'met' or 'slice' sampling should be used for the region-level effects 204 Lambda_configs : configuration options for the sampler for the region-level effects that will be used in the step method. 205 206 For options that can be in Lambda/Rho_configs, see: 207 spvcm.steps.Slice, spvcm.steps.Metropolis 208 truncation : dictionary 209 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: 210 Rho_min : minimum value allowed for response-level autoregressive coefficient 211 Rho_max : maximum value allowed for response-level autoregressive coefficient 212 Lambda_min : minimum value allowed for region-level autoregressive coefficient 213 Lambda_max : maximum value allowed for region-level autoregressive coefficient. 214 """ 215 def __init__(self, Y, X, M, Z=None, Delta=None, membership=None, 216 #data options 217 transform ='r', verbose=False, 218 n_samples=1000, n_jobs=1, 219 extra_traced_params = None, 220 priors=None, 221 configs=None, 222 starting_values=None, 223 truncation=None, 224 center=False, 225 scale=False, 226 tuning=0): 227 if X is None: 228 X = np.ones_like(Y) 229 center=False 230 scale=False 231 _, M = verify.weights(None, M, transform=transform) 232 self.M = M 233 Mmat = M.sparse 234 235 N,_ = X.shape 236 if Delta is not None: 237 J = Delta.shape[1] 238 elif membership is not None: 239 J = len(np.unique(membership)) 240 241 Delta, membership = verify.Delta_members(Delta, membership, N, J) 242 if Z is not None: 243 Z = Delta.dot(Z) 244 X = np.hstack((X,Z)) 245 if center: 246 X = verify.center(X) 247 if scale: 248 X = verify.scale(X) 249 250 X = verify.covariates(X) 251 252 self._verbose = verbose 253 254 super(Upper_SE, self).__init__(Y, X, Mmat, Delta, 255 n_samples=n_samples, 256 n_jobs = n_jobs, 257 extra_traced_params=extra_traced_params, 258 priors=priors, 259 configs=configs, 260 starting_values=starting_values, 261 truncation=truncation) 262