Skip to content

Commit 58e3964

Browse files
baddstatsbaddstats
committed
run Markov chain with sparse transition matrix
1 parent 2ec9bab commit 58e3964

12 files changed

Lines changed: 391 additions & 15 deletions

File tree

‎DESCRIPTION‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
Package: spatstat.sparse
2-
Version: 3.1-0
3-
Date: 2024-06-21
2+
Version: 3.1-0.001
3+
Date: 2026-04-10
44
Title: Sparse Three-Dimensional Arrays and Linear Algebra Utilities
55
Authors@R: c(person("Adrian", "Baddeley",
66
role = c("aut", "cre", "cph"),

‎NAMESPACE‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -53,6 +53,7 @@ export("print.sparse3Darray")
5353
export("quadform")
5454
export("rbindCompatibleDataFrames")
5555
export("representativeRows")
56+
export("runSparseMarkovChain")
5657
export("[.sparse3Darray")
5758
export("[<-.sparse3Darray")
5859
export("sparse3Darray")

‎NEWS‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,16 @@
11

2+
CHANGES IN spatstat.sparse VERSION 3.1-0.001
3+
4+
OVERVIEW
5+
6+
o Simulate a Markov chain with sparse transition matrix.
7+
8+
NEW FUNCTIONS
9+
10+
o runSparseMarkovChain
11+
Generate simulated realisations of a Markov chain
12+
using a sparse matrix representation of the transition probability matrix.
13+
214
CHANGES IN spatstat.sparse VERSION 3.1-0
315

416
OVERVIEW

‎R/sparse3Darray.R‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@
66
#' Copyright (c) Adrian Baddeley, Ege Rubak and Rolf Turner 2016-2020
77
#' GNU Public Licence >= 2.0
88
#'
9-
#' $Revision: 1.45 $ $Date: 2023/06/23 02:26:26 $
9+
#' $Revision: 1.46 $ $Date: 2026/01/21 06:26:39 $
1010
#'
1111

1212
sparse3Darray <- function(i=integer(0), j=integer(0), k=integer(0),
@@ -203,7 +203,7 @@ aperm.sparse3Darray <- function(a, perm=NULL, resize=TRUE, ...) {
203203
a$dim <- a$dim[perm]
204204
if(length(a$dimnames)==3) a$dimnames <- a$dimnames[perm]
205205
}
206-
class(a) <- c("sparse3Darray", class(a))
206+
class(a) <- unique(c("sparse3Darray", class(a)))
207207
return(a)
208208
}
209209

‎R/sparseMarkov.R‎

Lines changed: 95 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,95 @@
1+
#'
2+
#' sparseMarkov.R
3+
#'
4+
#' discrete-time finite-state Markov chain simulation
5+
#' using a sparse matrix representation of the transition matrix
6+
#'
7+
#' $Revision: 1.4 $ $Date: 2026/04/10 03:58:24 $
8+
#'
9+
#' Copyright (c) Adrian Baddeley 2026
10+
#' GNU Public Licence (>= 2.0)
11+
12+
13+
runSparseMarkovChain <- function(P, x0, nsteps, ...,
14+
result=c("history", "last"),
15+
check=TRUE,
16+
method=c("C", "interpreted")) {
17+
P <- as(P, "RsparseMatrix")
18+
if(!inherits(P, "dgRMatrix"))
19+
stop("Unable to convert P to a sparse matrix in row-major form",
20+
call.=FALSE)
21+
method <- match.arg(method)
22+
result <- match.arg(result)
23+
savehistory <- (result == "history")
24+
x0 <- as.integer(x0)
25+
nx <- length(x0)
26+
if(check) {
27+
ra <- range(P)
28+
if(ra[1L] < 0) stop("P contains negative entries", call.=FALSE)
29+
if(ra[2L] > 1) stop("P contains entries greater than 1", call.=FALSE)
30+
rs <- range(rowSums(P))
31+
if(max(abs(rs-1)) > sqrt(.Machine$double.eps))
32+
stop("Some of the rows of P do not sum to 1", call.=FALSE)
33+
rx <- range(x0)
34+
if(rx[1L] < 1 || rx[2L] > nrow(P))
35+
stop("Some indices in x0 are out of range", call.=FALSE)
36+
}
37+
## >>>>>>>>>>>> run chain <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
38+
switch(method,
39+
interpreted = {
40+
x <- x0
41+
if(savehistory) history <- x0
42+
for(istep in 1:nsteps) {
43+
xnew <- integer(0)
44+
## make one transition for each particle
45+
for(i in x) {
46+
p.i <- P[i, ]
47+
jj <- which(p.i > 0)
48+
p.ij <- p.i[jj]
49+
xnew <- c(xnew, sample(jj, 1L, prob=p.ij))
50+
}
51+
## current state of each particle
52+
x <- xnew
53+
## append to history (column = step, row = particle)
54+
if(savehistory) history <- cbind(history, x)
55+
}
56+
},
57+
C = {
58+
if(savehistory) {
59+
z <- .C(SP_rMCspMH,
60+
nrows = as.integer(dim(P)[1L]),
61+
nsparse = as.integer(length(P@x)),
62+
probval = as.double(P@x),
63+
colindex = as.integer(P@j),
64+
rowstart = as.integer(P@p),
65+
npoints = as.integer(nx),
66+
startpos = as.integer(x0 - 1L), # zero-based indices
67+
nsteps = as.integer(nsteps),
68+
history = as.integer(integer(nx * (nsteps + 1L))),
69+
PACKAGE = "spatstat.sparse")
70+
history <- matrix(z$history + 1L, nx, nsteps+1L, byrow=TRUE)
71+
} else {
72+
z <- .C(SP_rMCspMF,
73+
nrows = as.integer(dim(P)[1L]),
74+
nsparse = as.integer(length(P@x)),
75+
probval = as.double(P@x),
76+
colindex = as.integer(P@j),
77+
rowstart = as.integer(P@p),
78+
npoints = as.integer(nx),
79+
startpos = as.integer(x0 - 1L), # zero-based indices
80+
nsteps = as.integer(nsteps),
81+
endpos = as.integer(integer(nx)),
82+
PACKAGE = "spatstat.sparse")
83+
x <- z$endpos + 1L
84+
}
85+
})
86+
if(savehistory) {
87+
dimnames(history) <- NULL
88+
if(nx == 1) history <- history[1L, ]
89+
return(history)
90+
}
91+
return(x)
92+
}
93+
94+
95+

‎inst/doc/packagesizes.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,3 +19,4 @@ date version nhelpfiles nobjects ndatasets Rlines srclines
1919
"2023-03-12" "3.0-1" 15 48 0 2092 740
2020
"2023-10-24" "3.0-3" 15 48 0 2092 740
2121
"2024-06-21" "3.1-0" 15 48 0 2092 740
22+
"2026-04-10" "3.1-0.001" 16 49 0 2187 908

‎inst/info/packagesizes.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,3 +19,4 @@ date version nhelpfiles nobjects ndatasets Rlines srclines
1919
"2023-03-12" "3.0-1" 15 48 0 2092 740
2020
"2023-10-24" "3.0-3" 15 48 0 2092 740
2121
"2024-06-21" "3.1-0" 15 48 0 2092 740
22+
"2026-04-10" "3.1-0.001" 16 49 0 2187 908

‎man/macros/defns.Rd‎

Lines changed: 24 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -5,20 +5,21 @@
55
\newcommand{\ege}{Ege Rubak \email{rubak@math.aau.dk}}
66
\newcommand{\spatstatAuthors}{\adrian, \rolf and \ege}
77
\newcommand{\spatstatAuthorsComma}{\adrian, \rolf, \ege}
8-
%% Contributors with emails
9-
\newcommand{\pavel}{Pavel Grabarnik \email{pavel.grabar@issp.serpukhov.su}}
10-
\newcommand{\dominic}{Dominic Schuhmacher \email{dominic.schuhmacher@mathematik.uni-goettingen.de}, URL \code{http://dominic.schuhmacher.name/}}
11-
\newcommand{\wei}{Ang Qi Wei \email{aqw07398@hotmail.com}}
12-
\newcommand{\colette}{Marie-Colette van Lieshout \email{Marie-Colette.van.Lieshout@cwi.nl}}
13-
\newcommand{\rasmus}{Rasmus Plenge Waagepetersen \email{rw@math.auc.dk}}
8+
%% Contributors with emails (alphabetical order)
149
\newcommand{\abdollah}{Abdollah Jalilian \email{jalilian@razi.ac.ir}}
10+
\newcommand{\colette}{Marie-Colette van Lieshout \email{Marie-Colette.van.Lieshout@cwi.nl}}
11+
\newcommand{\dominic}{Dominic Schuhmacher \email{dominic.schuhmacher@mathematik.uni-goettingen.de}, URL \code{http://dominic.schuhmacher.name/}}
12+
\newcommand{\martinH}{Martin Hazelton \email{Martin.Hazelton@otago.ac.nz}}
13+
\newcommand{\mehdi}{Mehdi Moradi \email{m2.moradi@yahoo.com}}
1514
\newcommand{\ottmar}{Ottmar Cronie \email{ottmar@chalmers.se}}
15+
\newcommand{\pavel}{Pavel Grabarnik \email{pavel.grabar@issp.serpukhov.su}}
16+
\newcommand{\rasmus}{Rasmus Plenge Waagepetersen \email{rw@math.auc.dk}}
1617
\newcommand{\stephenEglen}{Stephen Eglen \email{S.J.Eglen@damtp.cam.ac.uk}}
17-
\newcommand{\mehdi}{Mehdi Moradi \email{m2.moradi@yahoo.com}}
18-
\newcommand{\yamei}{Ya-Mei Chang \email{yamei628@gmail.com}}
19-
\newcommand{\martinH}{Martin Hazelton \email{Martin.Hazelton@otago.ac.nz}}
18+
\newcommand{\suman}{Suman Rakshit \email{Suman.Rakshit@curtin.edu.au}}
2019
\newcommand{\tilman}{Tilman Davies \email{Tilman.Davies@otago.ac.nz}}
21-
% Names with accents
20+
\newcommand{\wei}{Ang Qi Wei \email{aqw07398@hotmail.com}}
21+
\newcommand{\yamei}{Ya-Mei Chang \email{yamei628@gmail.com}}
22+
% Names with accents (alphabetical order)
2223
\newcommand{\Bogsted}{\ifelse{latex}{\out{B\o gsted}}{Bogsted}}
2324
\newcommand{\Cramer}{\ifelse{latex}{\out{Cram\'er}}{Cramer}}
2425
\newcommand{\Francois}{\ifelse{latex}{\out{Fran\c{c}ois}}{Francois}}
@@ -36,13 +37,25 @@
3637
\newcommand{\Dominguez}{\ifelse{latex}{\out{Dom\'{\i}nguez}}{Dominguez}}
3738
\newcommand{\Rodriguez}{\ifelse{latex}{\out{Rodr\'{\i}guez}}{Rodriguez}}
3839
\newcommand{\Gonzalez}{\ifelse{latex}{\out{Gonz\'{a}lez}}{Gonzalez}}
40+
\newcommand{\Jimenez}{\ifelse{latex}{\out{Jim{\'e}nez}}{Jimenez}}
41+
\newcommand{\Birkhauser}{\ifelse{latex}{\out{Birkh{\"a}user}}{Birkhauser}}
3942
%% List of all Gibbs interactions
4043
\newcommand{\GibbsInteractionsList}{\code{\link[MPKG]{AreaInter}}, \code{\link[MPKG]{BadGey}}, \code{\link[MPKG]{Concom}}, \code{\link[MPKG]{DiggleGatesStibbard}}, \code{\link[MPKG]{DiggleGratton}}, \code{\link[MPKG]{Fiksel}}, \code{\link[MPKG]{Geyer}}, \code{\link[MPKG]{Hardcore}}, \code{\link[MPKG]{HierHard}}, \code{\link[MPKG]{HierStrauss}}, \code{\link[MPKG]{HierStraussHard}}, \code{\link[MPKG]{Hybrid}}, \code{\link[MPKG]{LennardJones}}, \code{\link[MPKG]{MultiHard}}, \code{\link[MPKG]{MultiStrauss}}, \code{\link[MPKG]{MultiStraussHard}}, \code{\link[MPKG]{OrdThresh}}, \code{\link[MPKG]{Ord}}, \code{\link[MPKG]{Pairwise}}, \code{\link[MPKG]{PairPiece}}, \code{\link[MPKG]{Penttinen}}, \code{\link[MPKG]{Poisson}}, \code{\link[MPKG]{Saturated}}, \code{\link[MPKG]{SatPiece}}, \code{\link[MPKG]{Softcore}}, \code{\link[MPKG]{Strauss}}, \code{\link[MPKG]{StraussHard}} and \code{\link[MPKG]{Triplets}}}
4144
%% List of interactions recognised by RMH code
4245
\newcommand{\rmhInteractionsList}{\code{\link[MPKG]{AreaInter}}, \code{\link[MPKG]{BadGey}}, \code{\link[MPKG]{DiggleGatesStibbard}}, \code{\link[MPKG]{DiggleGratton}}, \code{\link[MPKG]{Fiksel}}, \code{\link[MPKG]{Geyer}}, \code{\link[MPKG]{Hardcore}}, \code{\link[MPKG]{Hybrid}}, \code{\link[MPKG]{LennardJones}}, \code{\link[MPKG]{MultiStrauss}}, \code{\link[MPKG]{MultiStraussHard}}, \code{\link[MPKG]{PairPiece}}, \code{\link[MPKG]{Penttinen}}, \code{\link[MPKG]{Poisson}}, \code{\link[MPKG]{Softcore}}, \code{\link[MPKG]{Strauss}}, \code{\link[MPKG]{StraussHard}} and \code{\link[MPKG]{Triplets}}}
4346
%% Frequent references
4447
\newcommand{\baddrubaturnbook}{Baddeley, A., Rubak, E. and Turner, R. (2015) \emph{Spatial Point Patterns: Methodology and Applications with R}. Chapman and Hall/CRC Press. }
4548
%% Citations of recent articles that will change rapidly
46-
\newcommand{\baddchangclustersim}{Baddeley, A. and Chang, Y.-M. (2023) Robust algorithms for simulating cluster point processes. \emph{Journal of Statistical Computation and Simulation}. In Press. DOI \code{10.1080/00949655.2023.2166045}.}
49+
\newcommand{\baddchangclustersim}{Baddeley, A. and Chang, Y.-M. (2023) Robust algorithms for simulating cluster point processes. \emph{Journal of Statistical Computation and Simulation} \bold{93}, 1950--1975.}
50+
\newcommand{\smoothpcfpaper}{Baddeley, A., Davies, T.M. and Hazelton, M.L. (2025) An improved estimator of the pair correlation function of a spatial point process. \emph{Biometrika} \bold{111}, 2, article asaf021.}
51+
\newcommand{\smoothpcfpapercite}{Baddeley, Davies and Hazelton (2025)}
52+
%% ROC paper
53+
\newcommand{\rocketAuthors}{\adrian, \ege and \suman}
54+
\newcommand{\rocpapercite}{Baddeley et al (2025)}
55+
\newcommand{\rocpapercitealp}{Baddeley et al, 2025)}
56+
\newcommand{\rocpaper}{Baddeley, A., Rubak, E., Rakshit, S. and Nair, G. (2025) ROC curves for spatial point patterns and presence-absence data. \doi{10.48550/arXiv.2506.03414}.}
57+
%% pointweights
58+
\newcommand{\pointweightsAny}{A numeric or logical valued vector, matrix, data frame, pixel image, \code{function} or \code{expression}, or one of the strings \code{"x"} or \code{"y"} representing the cartesian coordinates. The \code{expression} may involve the variables \code{x,y,marks} representing the coordinates and marks of the point pattern}
59+
\newcommand{\pointweightsVector}{A numeric or logical valued vector, pixel image, \code{function} or \code{expression}, or one of the strings \code{"x"} or \code{"y"} representing the cartesian coordinates. The \code{expression} may involve the variables \code{x,y,marks} representing the coordinates and marks of the point pattern}
4760
4861

‎man/runSparseMarkovChain.Rd‎

Lines changed: 85 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,85 @@
1+
\name{runSparseMarkovChain}
2+
\alias{runSparseMarkovChain}
3+
\title{
4+
Run a Markov Chain with a Sparse Transition Matrix
5+
}
6+
\description{
7+
Generate a simulated realisation of a discrete-time finite-state
8+
Markov chain, using a sparse matrix representation of the transition
9+
probability matrix.
10+
}
11+
\usage{
12+
runSparseMarkovChain(P, x0, nsteps, \dots,
13+
result = c("history", "last"), check = TRUE,
14+
method = c("C", "interpreted"))
15+
}
16+
\arguments{
17+
\item{P}{
18+
Transition matrix. A sparse matrix supported by the
19+
\pkg{Matrix} package, or a \code{matrix} in the base \R package.
20+
}
21+
\item{x0}{
22+
Integer between \code{1} and \code{nrow(P)}
23+
specifying the initial state. Alternatively a vector of integers.
24+
}
25+
\item{nsteps}{
26+
Integer. Number of steps of the Markov chain.
27+
}
28+
\item{\dots}{
29+
Ignored.
30+
}
31+
\item{result}{
32+
Character string (partially matched) indicating whether to return
33+
the entire trajectory of the chain (\code{result="history"}, the
34+
default) or just the final state of the chain (\code{result="last"}).
35+
}
36+
\item{check}{
37+
Logical value specifying whether to check the validity of \code{P}
38+
and \code{x0}.
39+
}
40+
\item{method}{
41+
Character string (partially matched) indicating whether to use
42+
the \proglang{C} code implementation or an interpreted \R
43+
implementation. For testing purposes only.
44+
}
45+
}
46+
\details{
47+
If \code{x0} is a single integer, the Markov chain
48+
with transition probability matrix \code{P} will be
49+
run for \code{nsteps} steps starting from initial state \code{x0}.
50+
The result will be a vector of length \code{nsteps+1} giving the
51+
history of the chain after \code{0, 1, \dots, nsteps} steps
52+
(if \code{result="history"}, the default) or a
53+
single integer giving the final state of the chain
54+
(if \code{result="last"}).
55+
56+
If \code{x0} is a vector of integers, each entry of \code{x0}
57+
will be taken as the starting state for an independent copy of the
58+
Markov chain. Each copy will be run for \code{nsteps} steps.
59+
The result will be a matrix with \code{n=length(x0)} rows
60+
and \code{nsteps+1} columns giving the history of each chain
61+
(if \code{result="history"}, the default) or a
62+
vector of \code{n} integers giving the final state of each chain
63+
(if \code{result="last"}).
64+
65+
The matrix \code{P} should ideally be a sparse matrix
66+
in row major form (class \code{"dgRMatrix"}) but will be converted
67+
to this form.
68+
}
69+
\value{
70+
Integer vector or matrix.
71+
}
72+
\author{
73+
\adrian
74+
}
75+
\examples{
76+
M <- diag(5)
77+
M[1,3] <- 1
78+
M[2, 2:3] <- 1
79+
M[3, 2:4] <- 1
80+
M[4, 1:5] <- 1
81+
M[5, 1:5] <- 1
82+
P <- M/rowSums(M)
83+
runSparseMarkovChain(P, 2, 20)
84+
}
85+
\keyword{datagen}

‎src/init.c‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,8 @@ static const R_CMethodDef CEntries[] = {
2929
{"CwsumDsymouter", (DL_FUNC) &CwsumDsymouter, 5},
3030
{"Cwsumouter", (DL_FUNC) &Cwsumouter, 5},
3131
{"Cwsumsymouter", (DL_FUNC) &Cwsumsymouter, 5},
32+
{"rMCspMF", (DL_FUNC) &rMCspMF, 9},
33+
{"rMCspMH", (DL_FUNC) &rMCspMH, 9},
3234
{NULL, NULL, 0}
3335
};
3436

0 commit comments

Comments
 (0)