# How to represent mixture state with cgrad in plots.jl

Hi all,

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 = 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.

Please check if something like this might help:

``````plot(t, y, c=RGB.(eachcol(probs)...), lw=3, leg=false)
plot!(t, probs, c=[:red :green :blue])
``````

1 Like

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 = 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)
``````

The following code simply plots discretes colors, if that helps:

``````colors = [:red, :lime, :blue, :brown]
pal = palette(colors)
plot(t, y, line_z = argmax.(eachrow(probs)), c=pal, lw=3, leg=false)
plot!(t,probs, c=permutedims(colors))
``````

Adding a color-filled background may help too:

``````colors = [:red, :lime, :blue, :brown]
pal = palette(colors)
cprobs = -1 .+ 2*cumsum(probs, dims=2)
plot(t, -1 .+ zero(t), fillrange = cprobs[:,1], c=colors, fa=0.2, lw=0)
plot!(t, cprobs[:,1], fillrange = cprobs[:,2], c=colors, fa=0.2, lw=0)
plot!(t, cprobs[:,2], fillrange = cprobs[:,3], c=colors, fa=0.2, lw=0)
plot!(t, cprobs[:,3], fillrange = cprobs[:,4], c=colors, fa=0.3, lw=0)
plot!(t, y, line_z = argmax.(eachrow(probs)), c=pal, lw=4, leg=false, cb=false, widen=false)
``````

2 Likes

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.

I don’t know enough about RGB to know whether this makes sense, but it looks like I can make weighted sums of colors:

``````c1 = RGB(.5, .3, .2)
c2 = RGB(.7, .3,  .8)
c3 = .5 * c1 + .5 * c2
[c1,c2,c3]
``````

If this is a sensible way to mixture colors, I can use the state probabilities as weights. Do you think that makes sense?

Here is my best attempt:

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 = 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)

colors = RGB[]
push!(colors, RGB(.5, .2, .2))
push!(colors, RGB(.2, .3,  .8))
push!(colors, RGB(.3, .9,  .2))
push!(colors, RGB(.7, .5,  .1))

color_vec = probs * colors
plot(t, y, c=color_vec, lw=3, leg=false)
``````

This might work, but I am amenable to improvements.

1 Like

The only concern for 4 or more states is that we will be mixing colors that are not “linearly independent”?

1 Like

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,t[n], 4*n)
probs_fine = [itp(t,i) for t in tfine, i in ix]
``````
1 Like

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 = 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,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`.

t and y are shorter than probs_tfine.
Try: tfine and yfine=sin.(tfine)

Ok. Thank you. I see. In my case, t and y are empirical observations. I"ll see if I can interpolate t and y to make the vectors equal in length.

I try your code and since I am re-learning Statistics and Probability again this could be useful (Markov chain, hidden semi-Markov).

I want to ask on this code:

``````push!(colors, RGB(.2, .3,  .8))
``````

Is there any link to know the recipe to create a color based on this RGB pattern?

``````push!(colors, RGB(.5, .2, .2))
push!(colors, RGB(.2, .3,  .8))
push!(colors, RGB(.3, .9,  .2))
push!(colors, RGB(.7, .5,  .1))
``````

About 4 lines of the code above, is each one representing the first state, till the fourth state?

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])
``````