add randtrans and randperiodichmm

src/PeriodicHiddenMarkovModels.jl

@@ -29,7 +29,8 @@ export
     fit_mle,
     viterbi,
     # utilities.jl
-    n_to_t
+    n_to_t,
+    randPeriodicHMM
 
 include("utilities.jl")
 for fname in ["periodichmm.jl", "mle.jl", "likelihoods.jl"], foldername in ["periodic"]

src/utilities.jl

@@ -22,6 +22,7 @@ issquare(A::AbstractMatrix) = size(A, 1) == size(A, 2)
 """
     istransmat(A) -> Bool
 
+Function taken from [HMMBase.jl](https://maxmouchet.github.io/HMMBase.jl/stable/).
 Return true if `A` is square and its rows sums to 1.
 """
 istransmat(A::AbstractMatrix) =
@@ -73,4 +74,69 @@ function update_a!(a::AbstractVector, α::AbstractMatrix, β::AbstractMatrix)
     for i in OneTo(K)
         a[i] /= c
     end
-end
+end
+
+
+# TODO: should be in PeriodicHMM.jl or taken from HMM.jl
+
+"""
+    randtransmat([rng,] prior) -> Matrix{Float64}
+
+Function taken from [HMMBase.jl](https://maxmouchet.github.io/HMMBase.jl/stable/).
+Generate a transition matrix where each row is sampled from `prior`.  
+The prior must be a multivariate probability distribution, such as a
+Dirichlet distribution.
+
+**Arguments**
+- `prior::MultivariateDistribution`: distribution over the transition matrix rows.
+
+**Example**
+```julia
+A = randtransmat(Dirichlet([0.1, 0.1, 0.1]))
+```
+"""
+function randtransmat(rng::AbstractRNG, prior::MultivariateDistribution)
+    K = length(prior)
+    A = Matrix{Float64}(undef, K, K)
+    for i in OneTo(K)
+        A[i, :] = rand(rng, prior)
+    end
+    @check PeriodicHiddenMarkovModels.istransmat(A)
+    A
+end
+
+randtransmat(prior::MultivariateDistribution) = randtransmat(GLOBAL_RNG, prior)
+
+
+"""
+    randtransmat([rng, ]K, α = 1.0) -> Matrix{Float64}
+
+Function taken from [HMMBase.jl](https://maxmouchet.github.io/HMMBase.jl/stable/).
+Generate a transition matrix where each row is sampled from
+a Dirichlet distribution of dimension `K` and concentration
+parameter `α`.
+
+**Arguments**
+- `K::Integer`: number of states.
+- `α::Float64 = 1.0`: concentration parameter of the Dirichlet distribution.
+
+**Example**
+```julia
+A = randtransmat(4)
+```
+"""
+randtransmat(rng::AbstractRNG, K::Integer, α = 1.0) = randtransmat(rng, Dirichlet(K, α))
+
+randtransmat(K::Integer, args...) = randtransmat(GLOBAL_RNG, K, args...)
+
+function randPeriodicHMM(rng::AbstractRNG, K, T, D; ξ=ones(K) / K)
+    B_rand = Bernoulli.(rand(rng, K, T)) # completly random -> bad
+    Q_rand = zeros(K, K, T)
+    for t in 1:T
+        Q_rand[:, :, t] = randtransmat(rng, K) # completly random -> bad
+    end
+    hmm_random = PeriodicHMM(ξ, Q_rand, B_rand)
+    return hmm_random
+end
+
+randPeriodicHMM(K, T, D; ξ=ones(K) / K) = randPeriodicHMM(GLOBAL_RNG, K, T, D; ξ=ξ)