Using 2D Normal Mixtures for Spatial Point Patterns

This post describes the design and usage of the sppmix package, which defines classes and methods for spatial point pattern data and mixture models. You can install the package by this command.

devtools::install_github("wangyuchen/sppmix")

A two dimensional normal mixture is a basic building block of this package. We created a S3 class normmix and related methods for operations with normal mixture. In this post, you’ll learn everything that related to a normmix object. In first section, class normmix is introduced. When provided with intensity and a window, a normal mixture becomes an intensity surface. We defined the class intensity_surface and introduced it in the second section. In the end, we showed how to generate point pattern data from those two classes.

Normal Mixtures

Definition of a Normal Mixture

The class normmix is defined to be a 2D normal mixture, which consists the following parameters:

  1. Number of mixture components \(m\).
  2. Probability \(p\) for each component normal distribution.
  3. Mean vector \(\mu\) for each component.
  4. Covariance matrix \(\Sigma\) for each component.

Where \(m\) should be greater than zero and the sum of \(p\) should be \(1\).

Create a 2D Normal Mixture

There are two ways to create a normal mixture in sppmix. If you know all the parameters beforehand, you can use normmix() to create it. If you don’t know specifically what specific mixture you want, you can use rnormmix() to simulate a normal mixture to play with.

Generate Mixture from fixed parameter

In the first approach, you’ll need to know to pass ps, mus and sigmas to normmix(). Probability \(p\)’s should be passed in as a vector, whereas \(\mu\)’s and \(\Sigma\)’s in a list, each component of the list corresponding to the component in the mixture. The following example shows how to create a normal mixture with given parameter.

ps <- c(.3, .7)
mus <- list(c(0.2, 0.3), c(.8, .7))
sigmas <- list(.01 * diag(2), .01 * diag(2))
mix1 <- normmix(ps = ps, mus = mus, sigmas = sigmas)
mix1
## Normal Mixture with 2 component(s).

It will be of class normmix, and there’s some default S3 methods associated with it. If you want more details on the normal mixture, you can use the summary method. It will give you every piece of information about this normal mixture.

summary(mix1)
## Component 1 is centered at [ 0.2 , 0.3 ] with probability 0.3 
##  covariance matrix:
##        0.01    0.00
##        0.00    0.01
## Component 2 is centered at [ 0.8 , 0.7 ] with probability 0.7 
##  covariance matrix:
##        0.01    0.00
##        0.00    0.01

You can also plot it using a density plot.

plot(mix1)

Generate Mixture from Simulation

If you don’t have a specific normal mixture you want to work with (which is always the case), you may want to simulate a normal mixture. rnormmix() will allow you to simulate normal mixtures in different ways.

Generating a normal mixture requires generating covariance matrix \(\Sigma\)’s from a Wishart distribution. We use a Wishart distribution \(W_2(V, n)\) with scale matrix \(V = \sigma_0 I_2\) to generate \(\Sigma\)’s. \(\sigma_0\) is a scale parameter and \(n\) is the degree of freedom.

mix2 <- rnormmix(m = 3, sig0 = .1, df = 5)
summary(mix2)
## Component 1 is centered at [ 0.04 , 0.53 ] with probability 0.5744 
##  covariance matrix:
##      0.5061 -0.0364
##     -0.0364  1.0980
## Component 2 is centered at [ 0.52 , 0.98 ] with probability 0.3451 
##  covariance matrix:
##      0.2796 -0.3605
##     -0.3605  0.8875
## Component 3 is centered at [ 0.08 , 0.64 ] with probability 0.0804 
##  covariance matrix:
##      0.4273 -0.1017
##     -0.1017  0.6788

By default, normmix() will restrict the component mean within \((0, 1)\). You can change that by specify xlim and ylim parameters. There’s also another option that allows you to simulate not only parameters, but also number of components. When rand_m is true, it will generate a mixture with random number of component where the maximum possible component number is \(m\).

mix3 <- rnormmix(m = 5, sig0 = .1, df = 5, rand_m = TRUE, ylim = c(0, 5))
summary(mix3)
## Component 1 is centered at [ 0.39 , 3.7 ] with probability 1 
##  covariance matrix:
##      0.3131 -0.1716
##     -0.1716  0.4219

Intensity Surface

Create Intensity Surface

An intensity surface generated from the normal mixture is another data structure we used a lot in this package. The intensity_surface class inherits from the normmix class. Besides the parameters for normal mixture, it will need two more parameters, a region (window), and an average number of points over that region.

There’s also two ways to create an intensity surface. You can either specify the parameters directly or add those additional parameters to a existing normmix object.

Create Intensity Surface from Scratch

To create an intensity surface from scratch, you’ll call normmix() with two more additional parameters: the average number of points lambda and an window object of class spatstat::owin.

intsurf1 <- normmix(ps = c(.3, .7),
                    mus = list(c(0.2, 0.2), c(.8, .8)),
                    sigmas = list(.01*diag(2), .01*diag(2)),
                    lambda = 100,
                    win = square(1))
intsurf1
## Expected number of points over window: 100 
## window: rectangle = [0, 1] x [0, 1] units
## Normal Mixture with 2 component(s).

You may noticed that besides the information for normal mixture, it will print out lambda and window. You can also summary an intensity_surface object.

summary(intsurf1)
## Expected number of points over window: 100 
## window: rectangle = [0, 1] x [0, 1] units
## Normal Mixture with 2 component(s).
## Component 1 is centered at [ 0.2 , 0.2 ] with probability 0.3 
##  covariance matrix:
##        0.01    0.00
##        0.00    0.01
## Component 2 is centered at [ 0.8 , 0.8 ] with probability 0.7 
##  covariance matrix:
##        0.01    0.00
##        0.00    0.01

Create Intensity Surface from Normal Mixture

If you already have a normmix, for example a randomly generated normmix from rnormmix(), you may want to extend it to be an intensity_surface. Here you can use the to_int_surf() function to add those additional parameters to the normmix object.

to_int_surf(mix1, lambda = 100, win = square(1))
## Expected number of points over window: 100 
## window: rectangle = [0, 1] x [0, 1] units
## Normal Mixture with 2 component(s).

This function can also be used when you need to change lambda or window for an intensity_surface object.

to_int_surf(intsurf1, win = square(2))
## Expected number of points over window: 100 
## window: rectangle = [0, 2] x [0, 2] units
## Normal Mixture with 2 component(s).
to_int_surf(intsurf1, lambda = 50)
## Expected number of points over window: 50 
## window: rectangle = [0, 1] x [0, 1] units
## Normal Mixture with 2 component(s).

Generating Point Pattern from Normal Mixtures

All the normmix objects above are used to represent theoretical models, now we try to simulate some data from those normal mixtures. The function rsppmix() can simulate a point pattern from a normal mixture.

An intensity surface is the minimum requirement for generating a point pattern from normal mixture. It specifies the region and intensity, and it has the normal mixture density for us to sample from.

rsppmix() will return a point pattern in the class of c("sppmix", "ppp"). You may recognize that ppp is the class for point pattern in the spatstat package. We extended it to sppmix, which is our own class for a point pattern generated from normal mixture.

pp1 <- rsppmix(intsurf1)
pp1
## Planar point pattern: 98 points
## window: rectangle = [0, 1] x [0, 1] units

You may also use a normmix object and supply the other information to sppmix. In this package, you’ll see many functions that expect a intsurf object and have ... parameter. In most cases, you can provide a normmix object instead and provide lambda and win to ... parameter. This allows you to use a call-specific intensity or window instead of create an intensity surface every time.

rsppmix(mix1, lambda = 100, win = square(1))
## Planar point pattern: 93 points
## window: rectangle = [0, 1] x [0, 1] units
plot(pp1)

Now you have an intensity surface and the point generated from it. You can use plot to plot the intensity and also add points to it.

plot(intsurf1, pattern = pp1)

Or an interactive 3D surface plot using plot_mix_3d.

plot_mix_3d(intsurf1)