I am working with a model that could be classified as a hidden semi-Markov model. The hidden state is not observable, but it could be estimated from the model and displayed graphically. I am interested in representing the state as a color gradient, but I don’t know exactly how to do so. In the simplified code below, I have a sin function and use a markov chain to simulate transitions across 3 states. I would like to plot the sin function and color code the line as a gradient according to the state distribution.
using StatsBase, Plots
function markov(tmat, s, n)
obs = fill(0, n)
n_row = length(s)
r = sample(1:n_row, Weights(s))
obs[1] = r
for i in 2:n
r = sample(1:n_row, Weights(tmat[r,:]))
obs[i] = r
end
return obs
end
# transition matrix
tmat = [.99 .01 0.0;
0.0 .99 .01;
0.0 0.0 1.0]
# initial state vector
s = [1,0,0]
# time points
t = 0:.1:30
n = length(t)
# sin function at t
y = sin.(t)
# sample underlying states
samples = mapreduce(_ -> markov(tmat, s, n), hcat, 1:30)
# compute empirical estimate of state probabilities
probs = mapreduce(x -> mean(samples .== x, dims=2), hcat, 1:3)
plot(t, y, palette=cgrad([:blue,:green,:purple]), leg=false)
The line gradient might look somewhat similar to this example. For now, I have picked blue green and purple to correspond with states 1, 2 and 3. Perhaps there are better choices. Ultimately, I would like to combine this with a ribbon plot described here.
Very nice! Thank you again for your help. In the three state case, this is a very elegant use of state mixture probabilities and RGB.
I should have been more clear in my description that I need to generalize to more than three states (at least to four). I apologize for not making that clear. Do you know of a way to use cgrad or something similar to support four states?
I extended the example to four states below
Summary
using StatsBase, Plots
function markov(tmat, s, n)
obs = fill(0, n)
n_row = length(s)
r = sample(1:n_row, Weights(s))
obs[1] = r
for i in 2:n
r = sample(1:n_row, Weights(tmat[r,:]))
obs[i] = r
end
return obs
end
# transition matrix
tmat = [.99 .01 0.0 0.0;
0.0 .99 .01 0.0;
0.0 0.0 .99 .01;
0.0 0.0 0.0 1.0]
# initial state vector
s = [1,0,0,0]
# time points
t = 0:.1:30
n = length(t)
# sin function at t
y = sin.(t)
# sample underlying states
samples = mapreduce(_ -> markov(tmat, s, n), hcat, 1:30)
# compute empirical estimate of state probabilities
probs = mapreduce(x -> mean(samples .== x, dims=2), hcat, 1:4)
plot(t, y, c=RGB.(eachcol(probs)...), lw=3, leg=false)
Thanks for the ideas. I want to explore whether I can incorporate a gradient into your second example. If that does not work, I think the background fill you proposed would be a good alternative.
Yeah. That would indeed be a limitation. I believe this would mean that there could be cases where the mapping between state probability vector and color would not be unique.
@rafael.guerra, I was hoping to get your advice on a related issue. In the screenshot below, I use your first solution to add a color gradient. As you can see, the transitions are not always smooth. Do you know of any way to use SmoothingSplines with probs? Thanks!
For 2D cubic spline interpolation of probs “along” the time axis, you could try:
using Interpolations
ix = 1:size(probs,2)
itp = CubicSplineInterpolation((t, ix), probs)
tfine = LinRange(t[1],t[n], 4*n)
probs_fine = [itp(t,i) for t in tfine, i in ix]
Sorry if this is a stupid question, but I didn’t find an answering during my search. Is it possible to use interpolation when the length of probs_fin is larger than t? In the example below, I get an error.
Summary
using StatsBase, Plots
function markov(tmat, s, n)
obs = fill(0, n)
n_row = length(s)
r = sample(1:n_row, Weights(s))
obs[1] = r
for i in 2:n
r = sample(1:n_row, Weights(tmat[r,:]))
obs[i] = r
end
return obs
end
# transition matrix
tmat = [.99 .01 0.0;
0.0 .99 .01;
0.0 0.0 1.0]
# initial state vector
s = [1,0,0]
# time points
t = 0:.1:30
n = length(t)
# sin function at t
y = sin.(t)
# sample underlying states
samples = mapreduce(_ -> markov(tmat, s, n), hcat, 1:30)
# compute empirical estimate of state probabilities
probs = mapreduce(x -> mean(samples .== x, dims=2), hcat, 1:3)
plot(t, y, c=RGB.(eachcol(probs)...), lw=3, leg=false)
using Interpolations
ix = 1:size(probs,2)
itp = CubicSplineInterpolation((t, ix), probs)
tfine = LinRange(t[1],t[n], 4*n)
probs_fine = [itp(t,i) for t in tfine, i in ix]
plot(t, y, c=RGB.(eachcol(probs_fine)...), lw=3, leg=false)
Unfortunately, in my use case, I cannot increase the resolution of t.
I’m not aware of a general recipe in Plots to create a vector of colors. As Rafael pointed out, if you have more than three states, the vector of colors will not be linearly independent. As a consequence, there is not a unique mapping between color and underlying probability distribution. However, if you have three states, Rafael’s code below will have a unique mapping. His code is equivalent to colors = [RGB(1,0,0), RGB(0,1,0), RGB(0,0,1)]. This works out nicely because the implied basis is orthogonal and linearly independent.
plot(t, y, c=RGB.(eachcol(probs)...), lw=3, leg=false)
plot!(t, probs, c=[:red :green :blue])
