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