Add baseratematrix

src/matrices.jl

@@ -2,10 +2,14 @@
 # Licensed under the MIT License. See LICENSE in the project root.
 # ------------------------------------------------------------------
 
+# ----------
+# EMPIRICAL
+# ----------
+
 """
     ratematrix(data, var; [parameters])
 
-Computes the transition rate matrix for categorical variable
+Estimate the transition rate matrix for categorical variable
 `var` stored in geospatial `data` using a minimum lag `minlag`.
 
 See also [`probmatrix`](@ref), [`countmatrix`](@ref)
@@ -15,7 +19,7 @@ ratematrix(data::AbstractGeoTable, var; minlag=_defaultminlag(data)) = _ratematr
 """
     probmatrix(data, var; [parameters])
 
-Computes the transition probability matrix for categorical variable
+Estimate the transition probability matrix for categorical variable
 `var` stored in geospatial `data` using a minimum lag `minlag`.
 
 See also [`ratematrix`](@ref), [`countmatrix`](@ref)
@@ -25,7 +29,7 @@ probmatrix(data::AbstractGeoTable, var; minlag=_defaultminlag(data)) = _probmatr
 """
     countmatrix(data, var; [parameters])
 
-Computes the transition count matrix for categorical variable
+Estimate the transition count matrix for categorical variable
 `var` stored in geospatial `data` using a minimum lag `minlag`.
 
 See also [`ratematrix`](@ref), [`probmatrix`](@ref)
@@ -95,6 +99,46 @@ function _countmatrix(dom, vals, minlag)
   C
 end
 
+# ------------
+# THEORETICAL
+# ------------
+
+"""
+    baseratematrix(l, p)
+
+Transition rate matrix from mean lengths `l` and proportions `p`
+that assumes random transitions based on relative proportions.
+
+The transition rate for the pair `i -> j` is given by `-1 / l[i]`
+if `i == j` and by `(p[j] / (1 - p[i])) / l[i]` otherwise.
+
+## References
+
+* Carle et al 1998. [Conditional Simulation of Hydrofacies Architecture:
+  A Transition Probability/Markov Approach](https://doi.org/10.2110/sepmcheg.01.147)
+"""
+function baseratematrix(l, p)
+  nₗ = length(l)
+  nₚ = length(p)
+
+  # sanity checks
+  if nₗ != nₚ
+    throw(ArgumentError("lengths and proportions must have the same length"))
+  end
+  if !all(pᵢ -> 0 ≤ pᵢ ≤ 1, p)
+    throw(ArgumentError("proportions must be in interval [0, 1]"))
+  end
+
+  # Eq. 17 and 18 of Carle et al 1998.
+  map(Iterators.product(1:nₗ, 1:nₗ)) do (i, j)
+    if i == j
+      -1 / l[i]
+    else
+      (p[j] / (1 - p[i])) / l[i]
+    end
+  end
+end
+
 # -----------------
 # HELPER FUNCTIONS
 # -----------------

test/matrices.jl

@@ -31,4 +31,10 @@
   R = GeoStatsFunctions.ratematrix(gtb, "FACIES")
   @test size(R) == (15, 15)
   @test all(<(0 / u"m"), diag(R))
+
+  # base transition rate matrix
+  R = GeoStatsFunctions.baseratematrix([1,2,3]u"m", [0.2,0.5,0.3])
+  @test R == [         -1/1    0.5/(1-0.2)/1    0.3/(1-0.2)/1
+              0.2/(1-0.5)/2             -1/2    0.3/(1-0.5)/2
+              0.2/(1-0.3)/3    0.5/(1-0.3)/3             -1/3]u"m^-1"
 end