# How to simulate a random signal sequence change

**URL:** <https://discourse.julialang.org/t/how-to-simulate-a-random-signal-sequence-change/47964>\
**Category:** Statistics\
**Tags:** question\
**Created:** [October 8, 2020, 3:42am UTC](https://discourse.julialang.org/t/how-to-simulate-a-random-signal-sequence-change/47964 "2020-10-08T03:42:08Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![aachener](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aachener/32/18112_2.png) [@aachener](https://discourse.julialang.org/u/aachener)\
**Post date:** [October 8, 2020, 3:42am UTC](https://discourse.julialang.org/t/how-to-simulate-a-random-signal-sequence-change/47964/1 "2020-10-08T03:42:09Z")

</div>

I have got a simple problem to calculate Posterior after some signal disturbances.

Say the input signal is a sequence of 0 and 1 with 70% of chance is 1. And somehow we know that after the transmission, the output have 0 an 1 of equal chance. I want to calculate what the probability when we get 1 in output and it is truly 1 in the input.

I write the following code to implement the simulation. I want to randomly set seeds and see what happened to my simulation.

```julia
using Random, StatsBase
k = rand(1:2000);
Random.seed!(k);
t = sample(0:1,weights([0.3,0.7]),1000);
Random.seed!(k);
r = sample(0:1,weights([0.5,0.5]),1000);
Post = sum(t[r.==1])/sum(r[r.==1])
Prior =sum(r[t.==1])/sum(t[t.==1])

```

The result is not right, with Post always 1.0 and Prior around 0.7.

I am thinking of several potential problems:

1. I was conceptually wrong to implement the simulation this way. Maybe r is not the correct realization of t after disturbance, and I should modify t base on my knowledge of the disturbance instead of generating r directly.
2. I used sample() in a wrong way.
3. I calculated Post and Prior wrongly.

Please help me out on this simulation problem.

---

<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:** [October 8, 2020, 11:00pm UTC](https://discourse.julialang.org/t/how-to-simulate-a-random-signal-sequence-change/47964/2 "2020-10-08T23:00:08Z")

</div>

I am not an expert at all, but my guess would be:

```julia
Prior = sum(t[t.==1])/length(t)
Post = sum(t[r.==1 .& t.==1])/sum(t[t.==1])

```

and both should read ~0.7.  
I.e., the probability of the input being in state 1 given that the transmission was in state 1, should be 0.7 because the transmission is not changing the two-state probabilities. Seems analogous to the red-yellow cab problem well explained here [Bayes and the cabs](http://people.bu.edu/chamley/SL/SL-LN1.pdf).

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [October 9, 2020, 9:19am UTC](https://discourse.julialang.org/t/how-to-simulate-a-random-signal-sequence-change/47964/3 "2020-10-09T09:19:21Z")

</div>

I think the problem is not fully specified. You don’t link the output with the input, you only talk about their individual distributions:

So you have a probability vector before transmission `[0.7 0.3]`, and transmission procedure representable as 2x2 stochastic matrix `Q` and you tell us that `[0.7 0.3]*Q = [0.5 0.5]`. There are many `Q` which give you that.

```julia
julia> Q(λ) = [λ 1-λ; 5/3-7λ/3 1-( 5/3-7λ/3)]
Q (generic function with 1 method)

julia> Q(0.4)
2×2 Array{Float64,2}:
 0.4 0.6
 0.733333 0.266667

julia> [0.7 0.3]*Q(0.4)
1×2 Array{Float64,2}:
0.5 0.5

```

Here the entries of this particular `Q = [0.4 0.6; 0.733 0.267]` specify that a `1` gets corrupted to a 0 with probability 0.6 and a zero `0` gets corrupted to a `1` with probability `0.7333`.

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [October 9, 2020, 4:56pm UTC](https://discourse.julialang.org/t/how-to-simulate-a-random-signal-sequence-change/47964/4 "2020-10-09T16:56:21Z")

</div>

PS:

```julia
n = 10000
p = [0.7,0.3]
t = sample(1:-1:0, weights(p), n)
Q(λ) = [λ 1-λ; 5/3-7λ/3 1-( 5/3-7λ/3)]
λ = 0.4
r = [sample(1:-1:0, weights(vec(Q(λ)[i==1 ? 1 : 2, :]))) for i in t]

```

Then the conditional probability of being 1 given that that 1 was observed is computed and estimated as follows

```julia
julia> prob = (p[1]*Q(λ)[1,1])/(p[1]*Q(λ)[1,1] + p[2]*Q(λ)[2, 1])
0.56

julia> prob_est = dot(t,r)/n/mean(r)
0.5670391061452514

```

---

<div class="post-metadata">

**Author:** ![aachener](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aachener/32/18112_2.png) [@aachener](https://discourse.julialang.org/u/aachener)\
**Post date:** [October 10, 2020, 3:22am UTC](https://discourse.julialang.org/t/how-to-simulate-a-random-signal-sequence-change/47964/5 "2020-10-10T03:22:10Z")

</div>

Thank you, @mschauer, very clear explanation and helpful solution, too!

There are additional info on the disturbance process for us to model r.

It seems that independently generate r is not the correct solution. I will try as you suggested.