-
Notifications
You must be signed in to change notification settings - Fork 21
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
1 parent
1027739
commit d407124
Showing
4 changed files
with
1,356 additions
and
0 deletions.
There are no files selected for viewing
Large diffs are not rendered by default.
Oops, something went wrong.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,171 @@ | ||
{ | ||
"cells": [ | ||
{ | ||
"cell_type": "markdown", | ||
"metadata": {}, | ||
"source": [ | ||
"The ICAR (Intrinsic Conditional Autoregressive) model is a statistical model commonly used in spatial statistics to analyze and model spatial data. It is specifically used to account for spatial dependencies or autocorrelation in the data, where observations in close proximity tend to be more similar than those farther apart. The ICAR model is a type of conditional autoregressive model that assumes that the value of a location in a spatial dataset depends on the values of its neighboring locations.\n", | ||
"\n", | ||
"Here are the key components of the ICAR model:\n", | ||
"\n", | ||
"- Neighborhood Structure: The ICAR model requires a predefined neighborhood structure. This structure specifies which locations are considered neighbors of each other. Common choices include defining neighbors based on contiguity, distance thresholds, or some other criteria.\n", | ||
"\n", | ||
"- Conditional Autoregressive Structure: In the ICAR model, each location's value is assumed to be conditionally dependent on the values of its neighbors. This conditional dependency is often modeled as a weighted sum of the neighboring values, where the weights are determined by the neighborhood structure. The weights are typically constrained such that they sum to zero, ensuring stationarity of the model.\n", | ||
"\n", | ||
"- Intrinsic Condition: The \"intrinsic\" part of the ICAR model refers to the fact that the sum of the weights for each location is constrained to be zero. This constraint helps ensure identifiability of the model parameters.\n", | ||
"\n", | ||
"The ICAR model is useful in various spatial applications, such as spatial interpolation, spatial smoothing, and spatial prediction. It is commonly employed in Bayesian spatial modeling, where prior distributions are assigned to the model parameters, and the posterior distribution is obtained through Bayesian inference methods like Markov Chain Monte Carlo (MCMC) sampling.\n", | ||
"\n", | ||
"By accounting for spatial dependencies through the ICAR model, statisticians and researchers can better analyze and model spatial data while considering the inherent spatial correlation that often exists in such datasets." | ||
] | ||
}, | ||
{ | ||
"cell_type": "markdown", | ||
"metadata": {}, | ||
"source": [ | ||
"To implement the ICAR (Intrinsic Conditional Autoregressive) model in NumPyro, you can define a custom NumPyro model using the Pyro probabilistic programming framework, which NumPyro is built on. Here's an example of how you can write the ICAR model in NumPyro:" | ||
] | ||
}, | ||
{ | ||
"cell_type": "code", | ||
"execution_count": 4, | ||
"metadata": {}, | ||
"outputs": [ | ||
{ | ||
"name": "stderr", | ||
"output_type": "stream", | ||
"text": [ | ||
"sample: 100%|██████████| 2000/2000 [00:02<00:00, 779.69it/s, 3 steps of size 7.13e-08. acc. prob=0.91] \n" | ||
] | ||
} | ||
], | ||
"source": [ | ||
"import jax\n", | ||
"import numpy as np\n", | ||
"import numpyro\n", | ||
"import numpyro.distributions as dist\n", | ||
"import jax.numpy as jnp\n", | ||
"from numpyro.infer import MCMC, NUTS\n", | ||
"\n", | ||
"def icar_model(y, neighbors):\n", | ||
" n = len(y)\n", | ||
" alpha = numpyro.sample('alpha', dist.Normal(0, 1))\n", | ||
" sigma = numpyro.sample('sigma', dist.HalfCauchy(1))\n", | ||
" \n", | ||
" # Calculate the precision matrix\n", | ||
" tau = 1.0 / (sigma**2)\n", | ||
" precision_matrix = tau * (jnp.diag(jnp.sum(neighbors, axis=1)) - neighbors)\n", | ||
" \n", | ||
" # Define the conditional autoregressive structure\n", | ||
" with numpyro.plate('data', n):\n", | ||
" mu = alpha + jnp.matmul(precision_matrix, y)\n", | ||
" y = numpyro.sample('y', dist.MultivariateNormal(mu, precision_matrix), obs=y)\n", | ||
"\n", | ||
"# Example usage:\n", | ||
"# Define your spatial data 'y' and neighbor structure 'neighbors'\n", | ||
"y = np.array([1.2, 2.5, 3.8, 4.1, 5.2]) # Replace with your data\n", | ||
"neighbors = np.array([[0, 1, 0, 0, 0],\n", | ||
" [1, 0, 1, 0, 0],\n", | ||
" [0, 1, 0, 1, 0],\n", | ||
" [0, 0, 1, 0, 1],\n", | ||
" [0, 0, 0, 1, 0]]) # Replace with your neighbor structure\n", | ||
"\n", | ||
"# Run MCMC to infer the parameters\n", | ||
"nuts_kernel = NUTS(icar_model)\n", | ||
"mcmc = MCMC(nuts_kernel, num_samples=1000, num_warmup=1000)\n", | ||
"mcmc.run(rng_key=jax.random.PRNGKey(0), y=y, neighbors=neighbors)\n", | ||
"\n", | ||
"# Extract posterior samples\n", | ||
"samples = mcmc.get_samples()\n" | ||
] | ||
}, | ||
{ | ||
"cell_type": "code", | ||
"execution_count": 10, | ||
"metadata": {}, | ||
"outputs": [ | ||
{ | ||
"data": { | ||
"text/plain": [ | ||
"(1000, 2)" | ||
] | ||
}, | ||
"execution_count": 10, | ||
"metadata": {}, | ||
"output_type": "execute_result" | ||
} | ||
], | ||
"source": [ | ||
"import numpyro.distributions as dist\n", | ||
"import jax\n", | ||
"import jax.numpy as jnp\n", | ||
"from numpyro.infer import MCMC, NUTS\n", | ||
"\n", | ||
"# Define the mean and covariance matrix for the multivariate normal distribution\n", | ||
"mean = jnp.array([0.0, 0.0]) # Mean vector\n", | ||
"cov_matrix = jnp.array([[1.0, 0.5], [0.5, 2.0]]) # Covariance matrix\n", | ||
"\n", | ||
"# Create a MultivariateNormal distribution\n", | ||
"mvn = dist.MultivariateNormal(mean, cov_matrix)\n", | ||
"\n", | ||
"# Sample from the distribution\n", | ||
"samples = mvn.sample(key=jax.random.PRNGKey(0), sample_shape=(1000,))\n", | ||
"\n", | ||
"samples.shape\n", | ||
"\n", | ||
"# Now, 'samples' contains 1000 samples from the multivariate normal distribution\n" | ||
] | ||
}, | ||
{ | ||
"cell_type": "markdown", | ||
"metadata": {}, | ||
"source": [ | ||
"In this code:\n", | ||
"\n", | ||
"- We define the ICAR model as a Python function icar_model that takes y (the spatial data) and neighbors (the neighbor structure) as arguments.\n", | ||
"\n", | ||
"- Inside the model, we sample the parameters alpha (the intercept) and sigma (the spatial correlation parameter).\n", | ||
"\n", | ||
"- We use a plate to handle multiple data points if needed. The core of the model is the calculation of mu, which represents the conditional expectation based on the ICAR structure.\n", | ||
"\n", | ||
"- Finally, we sample y from a Normal distribution with mean mu and observation noise of 1.\n", | ||
"\n", | ||
"- We provide example usage of the model, where you should replace the y and neighbors arrays with your actual data and neighbor structure.\n", | ||
"\n", | ||
"- We run MCMC using the NUTS sampler to infer the posterior distribution of the model parameters.\n", | ||
"\n", | ||
"- After running MCMC, you can extract posterior samples from the mcmc object.\n", | ||
"\n", | ||
"Make sure to adapt this code to your specific dataset and spatial structure as needed." | ||
] | ||
}, | ||
{ | ||
"cell_type": "markdown", | ||
"metadata": {}, | ||
"source": [] | ||
} | ||
], | ||
"metadata": { | ||
"colab": { | ||
"provenance": [] | ||
}, | ||
"kernelspec": { | ||
"display_name": "Python 3", | ||
"name": "python3" | ||
}, | ||
"language_info": { | ||
"codemirror_mode": { | ||
"name": "ipython", | ||
"version": 3 | ||
}, | ||
"file_extension": ".py", | ||
"mimetype": "text/x-python", | ||
"name": "python", | ||
"nbconvert_exporter": "python", | ||
"pygments_lexer": "ipython3", | ||
"version": "3.9.7" | ||
} | ||
}, | ||
"nbformat": 4, | ||
"nbformat_minor": 0 | ||
} |
Large diffs are not rendered by default.
Oops, something went wrong.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,4 @@ | ||
jupyter-book | ||
matplotlib | ||
numpy | ||
ghp-import |