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

**URL:** <https://discourse.julialang.org/t/how-to-represent-mixture-state-with-cgrad-in-plots-jl/84706>\
**Category:** General Usage\
**Tags:** question, plotting, plots\
**Created:** [July 24, 2022, 12:23pm UTC](https://discourse.julialang.org/t/how-to-represent-mixture-state-with-cgrad-in-plots-jl/84706 "2022-07-24T12:23:33Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [July 24, 2022, 12:23pm UTC](https://discourse.julialang.org/t/how-to-represent-mixture-state-with-cgrad-in-plots-jl/84706/1 "2022-07-24T12:23:33Z")

</div>

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.

```julia
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](https://user-images.githubusercontent.com/5888590/103132482-c69af100-46b5-11eb-9e6e-0c642a87ea49.png). 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](https://discourse.julialang.org/t/add-gradient-to-ribbon-based-on-density/62081).

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [July 24, 2022, 9:44pm UTC](https://discourse.julialang.org/t/how-to-represent-mixture-state-with-cgrad-in-plots-jl/84706/2 "2022-07-24T21:44:23Z")

</div>

Please check if something like this might help:

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

```

 ![Plots_gr_RGB_line_color](https://global.discourse-cdn.com/julialang/original/3X/7/9/7941d8e95df8442a3d795aa6e561ecc70f05550f.png)

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [July 24, 2022, 10:51pm UTC](https://discourse.julialang.org/t/how-to-represent-mixture-state-with-cgrad-in-plots-jl/84706/3 "2022-07-24T22:51:13Z")

</div>

> [@rafael.guerra](#):
>
> `c=RGB.(eachcol(probs)...)`

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**
>
> ```julia
> 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)
> 
> ```

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [July 24, 2022, 11:50pm UTC](https://discourse.julialang.org/t/how-to-represent-mixture-state-with-cgrad-in-plots-jl/84706/4 "2022-07-24T23:50:19Z")

</div>

> [@Christopher\_Fisher](#):
>
> to support four states

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

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

```

 ![Plots_gr_line_color_4states](https://global.discourse-cdn.com/julialang/original/3X/1/6/160a0610cc9df5d6395c19665b833e4b32e3ce75.png)

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [July 25, 2022, 12:16am UTC](https://discourse.julialang.org/t/how-to-represent-mixture-state-with-cgrad-in-plots-jl/84706/5 "2022-07-25T00:16:01Z")

</div>

Adding a color-filled background may help too:

```julia
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[1], fa=0.2, lw=0)
plot!(t, cprobs[:,1], fillrange = cprobs[:,2], c=colors[2], fa=0.2, lw=0)
plot!(t, cprobs[:,2], fillrange = cprobs[:,3], c=colors[3], fa=0.2, lw=0)
plot!(t, cprobs[:,3], fillrange = cprobs[:,4], c=colors[4], fa=0.3, lw=0)
plot!(t, y, line_z = argmax.(eachrow(probs)), c=pal, lw=4, leg=false, cb=false, widen=false)

```

 ![Plots_gr_line_color_4states_filled_background](https://global.discourse-cdn.com/julialang/original/3X/1/c/1c3c56846b0f4c9c9733abca059a8aa3f0944631.png)

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [July 25, 2022, 1:01am UTC](https://discourse.julialang.org/t/how-to-represent-mixture-state-with-cgrad-in-plots-jl/84706/6 "2022-07-25T01:01:27Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [July 25, 2022, 1:15am UTC](https://discourse.julialang.org/t/how-to-represent-mixture-state-with-cgrad-in-plots-jl/84706/7 "2022-07-25T01:15:33Z")

</div>

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

```julia
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?

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [July 25, 2022, 1:26am UTC](https://discourse.julialang.org/t/how-to-represent-mixture-state-with-cgrad-in-plots-jl/84706/8 "2022-07-25T01:26:54Z")

</div>

Here is my best attempt:

> **Summary**
>
> ```julia
> 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)
> 
> 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.

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [July 25, 2022, 9:55am UTC](https://discourse.julialang.org/t/how-to-represent-mixture-state-with-cgrad-in-plots-jl/84706/9 "2022-07-25T09:55:37Z")

</div>

> [@Christopher\_Fisher](#):
>
> Do you think that makes sense?

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

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [July 25, 2022, 10:05am UTC](https://discourse.julialang.org/t/how-to-represent-mixture-state-with-cgrad-in-plots-jl/84706/10 "2022-07-25T10:05:06Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [July 28, 2022, 4:10pm UTC](https://discourse.julialang.org/t/how-to-represent-mixture-state-with-cgrad-in-plots-jl/84706/11 "2022-07-28T16:10:13Z")

</div>

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

@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!

 ![gradient](https://global.discourse-cdn.com/julialang/original/3X/f/5/f5b6ead3fbb35f1e2c2c495e3bfde44f50c1617e.png)

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [July 28, 2022, 9:04pm UTC](https://discourse.julialang.org/t/how-to-represent-mixture-state-with-cgrad-in-plots-jl/84706/12 "2022-07-28T21:04:36Z")

</div>

For 2D cubic spline interpolation of `probs` “along” the time axis, you could try:

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

```

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [July 28, 2022, 11:54pm UTC](https://discourse.julialang.org/t/how-to-represent-mixture-state-with-cgrad-in-plots-jl/84706/13 "2022-07-28T23:54:15Z")

</div>

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**
>
> ```julia
> 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`.

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [July 28, 2022, 11:58pm UTC](https://discourse.julialang.org/t/how-to-represent-mixture-state-with-cgrad-in-plots-jl/84706/14 "2022-07-28T23:58:46Z")

</div>

> [@Christopher\_Fisher](#):
>
> `plot(t, y, c=RGB.(eachcol(probs_fine)...), lw=3, leg=false)`

t and y are shorter than probs\_tfine.  
Try: tfine and yfine=sin.(tfine)

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [July 29, 2022, 12:14am UTC](https://discourse.julialang.org/t/how-to-represent-mixture-state-with-cgrad-in-plots-jl/84706/15 "2022-07-29T00:14:38Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Freya\_the\_Goddess](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/freya_the_goddess/32/36835_2.png) [@Freya\_the\_Goddess](https://discourse.julialang.org/u/Freya_the_Goddess)\
**Post date:** [August 2, 2022, 4:40am UTC](https://discourse.julialang.org/t/how-to-represent-mixture-state-with-cgrad-in-plots-jl/84706/16 "2022-08-02T04:40:43Z")

</div>

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:

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

```

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

```julia
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?

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [August 3, 2022, 11:21am UTC](https://discourse.julialang.org/t/how-to-represent-mixture-state-with-cgrad-in-plots-jl/84706/17 "2022-08-03T11:21:32Z")

</div>

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.

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

```
