1
$\begingroup$

I've simulated an N-by-N raster in the following way:

  1. define a set $S$ containing a finite number $|S| = K$ of possible raster values (in my simulation, $K=3$ and the elements of $S$ are land uses / land covers)
  2. Initialize an empty $N$-by-$N$ array of raster values; draw the $1, 1$ entry of the raster uniformly at random from $S$. I'll refer to the raster entries as $X_{i,j}$. At this point we have drawn a value for $X_{1, 1}$
  3. Populate the first row of the raster using a conditional distribution $\Pr[X_{i,j} \,|\, X_{i,j-1}]$, i.e. populate $X_{1, j}$ for $j=2, 3, \ldots, N$
  4. Use the same conditional distribution as above to populate the first column of the raster (i.e. $\Pr[X_{i,j} \,|\, X_{i-1,j}]$ is the same distribution as $\Pr[X_{i,j} \,|\, X_{i,j-1}]$); this allows us to populate $X_{i, 1}$ for $i=2, 3, \ldots, N$
  5. Populate the remainder of the raster using a distribution that conditions on two neighbors: $\Pr[X_{i,j} \,|\, X_{i,j-1}, X_{i-1,j}]$. I assume that this distribution satisfies a symmetry property: the probability does not change if you switch the values of the two neighbors, i.e. $\Pr[X_{i,j} \,|\, X_{i,j-1} = r, X_{i-1,j} = s]$ = $\Pr[X_{i,j} \,|\, X_{i,j-1} = s, X_{i-1,j} = r]$, for any $s, r \in S$.

Those steps generate an $N$-by-$N$ raster whose cells take values in $S$.

I have several questions:

  • The primitives of the model are the distributions $\Pr[X_{i,j} \,|\, X_{i,j-1}]$ and $\Pr[X_{i,j} \,|\, X_{i,j-1}, X_{i-1,j}]$. Do these objects imply any sort of stationary distribution over raster values $X_{i, j}$, and over the joint distribution of $X_{i, j}$ and its neighbors to the north, south, east, and west, i.e. the joint distribution of $\left(X_{i, j}, X_{i+1, j}, X_{i, j+1}, X_{i-1, j}, X_{i, j-1}\right)$?
  • If the primitives do imply a stationary distribution, can the stationary distribution be calculated using the primitives (the same way we can calculate the stationary distribution of an one-dimensional, first-order Markov chain using an eigenvector decomposition of its transition matrix)?
  • If there is a stationary joint distribution of $X_{i, j}$ and its neighbors to the north, south, east, and west, is that distribution "invariant to rotation", by which I mean that rotations of the neighbors around the central pixel do not change the probability? I suspect the answer is no: see the simulated raster below, which shows a sort of "streak" pattern from the 1, 1 pixel towards the top right (i.e., streaks moving in the direction in which I populated the raster in step 5).

Are there any references that answer these questions? I think what I'm simulating is a Markov random field on a grid, but I'm not sure what the proper terminology is.

R simulation:

library(ggplot2)

## Spatial dependence: distribution of ij raster pixel depends on its neighbors to north and west

S_sorted <- sort(c("crops", "forest", "pasture"))  # Set of land cover states
stopifnot(length(S_sorted) == 3)  # Assumed throughout

## Matrix with ij entry Pr[ S = s_j | S^{west} = s_i ]
## Symmetry assumption: this matrix is also used for Pr[ S = s_j | S^{north} = s_i ]
conditional_one_neighbor <- rbind(c(0.80, 0.10, 0.10),
                                  c(0.05, 0.90, 0.05),
                                  c(0.20, 0.10, 0.70))
stopifnot(isTRUE(all.equal(rowSums(conditional_one_neighbor),
                           rep(1.0, nrow(conditional_one_neighbor)))))  # Valid conditional probabilities

## Array with ijk entry Pr[ S = s_k | S^{north} = s_i, S^{west} = s_j ]
## Symmetry assumption: ijk entry also gives Pr[ S = s_k | S^{north} = s_j, S^{west} = s_i ],
## i.e. ijk entry and jik entry are always identical
conditional_two_neighbors <- array(0, dim=c(3, 3, 3))
conditional_two_neighbors[, , 1] <- rbind(c(0.95, 0.49, 0.29),
                                          c(0.49, 0.01, 0.01),
                                          c(0.29, 0.01, 0.02))  # Probability of crops
conditional_two_neighbors[, , 2] <- rbind(c(0.04, 0.50, 0.01),
                                          c(0.50, 0.95, 0.40),
                                          c(0.01, 0.40, 0.08))  # Probability of forest
conditional_two_neighbors[, , 3] <- rbind(c(0.01, 0.01, 0.70),
                                          c(0.01, 0.04, 0.59),
                                          c(0.70, 0.59, 0.90))  # Probability of pasture
for(k in seq_along(S_sorted)) {
    stopifnot(all(conditional_two_neighbors[, , k] == t(conditional_two_neighbors[, , k])))  # Symmetry assumption
}

## Sum over ij should be 1, i.e. valid probability distribution
stopifnot(isTRUE(all.equal(apply(conditional_two_neighbors, MARGIN=c(1, 2), FUN=sum),
                           matrix(1, dim(conditional_two_neighbors)[1], dim(conditional_two_neighbors)[2]))))

nrow_raster <- 500
raster <- matrix(NA, nrow_raster, nrow_raster)  # Simulate a square land cover raster
raster[1, 1] <- sample(S_sorted, size=1)  # First pixel chosen uniformly at random

## Populate (rest of) first row using conditional_one_neighbor
for(i in seq(2, nrow_raster)) {
    neighbor_index <- which(S_sorted == raster[1, i-1])  # Index into S_sorted
    raster[1, i] <- sample(S_sorted, size=1, prob=conditional_one_neighbor[neighbor_index, ])
}

## Populate (rest of) first column using conditional_one_neighbor
for(i in seq(2, nrow_raster)) {
    neighbor_index <- which(S_sorted == raster[i-1, 1])  # Index into S_sorted
    raster[i, 1] <- sample(S_sorted, size=1, prob=conditional_one_neighbor[neighbor_index, ])
}

## Populate rest of raster using conditional_two_neighbors (conditional on north and west neighbors)
for(i in seq(2, nrow_raster)) {
    for(j in seq(2, nrow_raster)) {
        neighbor_north_index <- which(S_sorted == raster[i-1, j])
        neighbor_west_index <- which(S_sorted == raster[i, j-1])
        raster[i, j] <- sample(S_sorted,
                               size=1,
                               prob=conditional_two_neighbors[neighbor_north_index, neighbor_west_index, ])
    }
}

raster_df <- expand.grid(x=seq_len(nrow_raster),
                         y=seq_len(nrow_raster))
raster_df$land_cover=as.vector(raster)
dim(raster_df)

## Raster origin (1, 1) is at bottom-left
p <- (ggplot(raster_df, aes(x=x, y=y, fill=land_cover)) +
      geom_raster() +
      scale_fill_manual("land cover", values=c("crops"="#ff7f00",
                                               "forest"="#4daf4a",
                                               "pasture"="#a6cee3")))
p
ggsave("simulated_raster.png", width=12, height=10)

## Convergence to stationary distribution as we move away from the origin, to the top-right?

table(raster_df$land_cover) / nrow(raster_df)  # Empirical marginal distribution

## Distribution conditional on one neighbor (empirical as above)
tables_list <- lapply(seq_len(nrow_raster), function(i) {
    return(table(raster[i, seq(1, nrow_raster-1)], raster[i, seq(2, nrow_raster)]))
})
tables_sum <- Reduce("+", tables_list)
joint_distrib <- tables_sum / sum(tables_sum)
empirical_conditional_one_neighbor <- joint_distrib / matrix(rep(rowSums(joint_distrib), 3), 3, 3)
rowSums(empirical_conditional_one_neighbor)  # Sanity check, all 1

## Compare (empirical) empirical_conditional_one_neighbor to conditional_one_neighbor: they're clearly different,
## one is an initial distribution and the other is something like a stationary distribution
round(empirical_conditional_one_neighbor, 3)
round(conditional_one_neighbor, 3)

simulated raster

When I say that I see "streaks" in the figure above, what I mean is that there are many clumps of forest that are long in the [down left]-to-[up right] direction, and not many clumps of forest that are long in the [up left]-to-[down right] direction. My intuition is that this is a consequence of the way that I simulated the raster, and that it means that, if there is a stationary joint distribution of $\left(X_{i, j}, X_{i+1, j}, X_{i, j+1}, X_{i-1, j}, X_{i, j-1}\right)$, i.e. of a pixel and its four neighbors, that distribution is not invariant to rotation. Is that correct?

Edit: has anyone reading this question read Introduction to Markov Random Fields by Andrew Blake and Pushmeet Kohli (chapter 1 available at https://mitpress.mit.edu/sites/default/files/titles/content/9780262015776_sch_0001.pdf)? If so, do you know whether I can find an answer to my question(s) in that book?

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.