# Help with Turing.jl model using Categorical, Normal and ordered constraint

**URL:** <https://discourse.julialang.org/t/help-with-turing-jl-model-using-categorical-normal-and-ordered-constraint/64765>\
**Category:** Probabilistic Programming\
**Tags:** turing\
**Created:** [July 16, 2021, 6:09pm UTC](https://discourse.julialang.org/t/help-with-turing-jl-model-using-categorical-normal-and-ordered-constraint/64765 "2021-07-16T18:09:49Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Hakan\_Kjellerstrand](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hakan_kjellerstrand/32/8447_2.png) [@Hakan\_Kjellerstrand](https://discourse.julialang.org/u/Hakan_Kjellerstrand)\
**Post date:** [July 16, 2021, 6:09pm UTC](https://discourse.julialang.org/t/help-with-turing-jl-model-using-categorical-normal-and-ordered-constraint/64765/1 "2021-07-16T18:09:50Z")

</div>

I’m trying to model AllenDowney’s three problems “How tall is A?” ([How tall is A? – Probably Overthinking It](https://www.allendowney.com/blog/2018/10/18/how-tall-is-a/) ):

```julia
1) Suppose you meet an adult resident of the U.S. who is 170 cm tall.
     What is the probability that they are male?

2) Suppose I choose two U.S. residents at random and A is taller than B. 
     How tall is A?

3) In a room of 10 randomly chosen U.S. residents, A is the second tallest. 
    How tall is A? 
    And what is the probability that A is male?*

```

The first two problems was easy to model in Turing.jl ([http://hakank.org/julia/turing/how\_tall\_is\_a.jl](http://hakank.org/julia/turing/how_tall_is_a.jl) ). Here’s the model of the second problem:

```julia
using Turing

@model function how_tall_is_2()
    male = 1
    female = 2
    gender = tzeros(2)
    height = tzeros(2)
    for p in 1:2
        gender[p] ~ Categorical(simplex([0.49,0.51]))
        height[p] ~ gender[p] == male ? Normal(178,7.7) : Normal(163,7.3)
    end
    true ~ Dirac(height[1] > height[2])    
end

model = how_tall_is_2()
chains = sample(model, PG(15), 1_000)
display(chains)

```

And it give the expected heights of:

```julia
   ...
   height[1] 176.4384 9.0135 0.2850 0.1777 824.6929 1.0007 179.5152
   height[2] 164.5469 8.4415 0.2669 0.2017 828.7247 0.9998 180.3928
   ...

```

However, the **third problem** , were we have 10 people instead of 2, is much harder to get correct. Here’s my model (about the same as the above, but `n` is now 10):

```julia
using Turing

@model function how_tall_is_3()
    male = 1
    female = 2
    n = 2
    gender = tzeros(n)
    height = tzeros(n)
    for p in 1:n
        gender[p] ~ Categorical(simplex([0.49,0.51]))
        height[p] ~ gender[p] == male ? Normal(178,7.7) : Normal(163,7.3)
    end

    for p in 1:n-1
        true ~ Dirac(height[p] > height[p+1])
        # height[p] > height[p+1] || begin Turing.@addlogprob! -Inf; end        
    end
 
end

model = how_tall_is_3()
chains = sample(model, MH(), 10_000)
# chains = sample(model, PG(15), 10_000)
# chains = sample(model, SMC(), 10_000)
# chains = sample(model, IS(), 10_000)
display(chains)

```

The correct answer should be that A (the second tallest person) is about 181.61 cm (and the tallest person about 186 cm).

But this model don’t give any ordered results. Here’s an (representative) example using MH() where the heights are all over the places instead of nicely ordered.

```julia
   height[1] 176.0685  
   height[2] 184.1998  
   height[3] 165.4559  
   height[4] 161.9992  
   height[5] 182.4204  
   height[6] 175.6319  
   height[7] 182.7167  
   height[8] 161.5589  
   height[9] 183.9556  
   height[10] 165.1349 
   ...

```

I’ve tested different samplers, e.g. HM(), PG(),SMC(), and IS() but all give some wrong answer, either with different but unordered heights, or heights that don’t differ much, e.g. PG() give all heights about 170 cm (i,e, the same as using Prior()) .

I also tested to use Gibbs mixing sampler, e.g. `Gibbs(MH(:gender),NUTS(1000,0.65,:height))` but it throws errors; probably because the parameters are not `:gender` and `:height`, but arrays of genders and heights and I’m not sure how to state this.

Also, as shown in the model, I tested both using `Dirac` and `Turing.@addlogprob!` for observing the order, but there’s no difference between these two.

I guess that the real problem is the mix of `Categorical` and `Normal` and that these samplers cannot handle the “observation” of ordering properly.

Is there a (preferably simple) fix for this, or perhaps the problem should be modeled in a completely different way?

---

<div class="post-metadata">

**Author:** ![sethaxen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sethaxen/32/35604_2.png) [@sethaxen](https://discourse.julialang.org/u/sethaxen)\
**Post date:** [July 16, 2021, 8:28pm UTC](https://discourse.julialang.org/t/help-with-turing-jl-model-using-categorical-normal-and-ordered-constraint/64765/2 "2021-07-16T20:28:37Z")

</div>

Here’s a model that does what you want:

```julia
using Turing

@model function how_tall_is_3(
    n,
    pfemale,
    h_mean_female,
    h_mean_male,
    h_std_female,
    h_std_male,
    height_A,
    isfemale_A,
    ::Type{T} = Float64,
) where {T}
    isfemale ~ filldist(Bernoulli(pfemale), n)

    height = Vector{T}(undef, n)
    for i in 1:n
        if isfemale[i]
            height[i] ~ Normal(h_mean_female, h_std_female)
        else
            height[i] ~ Normal(h_mean_male, h_std_male)
        end
    end
    id_A = sortperm(height)[n-1]
    height_A ~ Dirac(height[id_A])
    isfemale_A ~ Dirac(isfemale[id_A])
end

# set "observations" of A to `missing`, so that we generate them
mod = how_tall_is_3(10, 0.51, 163, 178, 7.3, 7.7, missing, missing)
# no inference necessary! We just sample directly from the prior-predictive distribution
chn = sample(mod, Prior(), 10_000)

```

Output:

```julia
Chains MCMC chain (10000×23×1 Array{Float64, 3}):

Iterations = 1:1:10000
Number of chains = 1
Samples per chain = 10000
Wall duration = 0.66 seconds
Compute duration = 0.66 seconds
parameters = height[10], height[1], height[2], isfemale[10], isfemale[5], isfemale[7], isfemale[4], isfemale[6], height_A, isfemale[8], height[7], height[8], isfemale[9], isfemale[3], height[4], height[6], height[9], isfemale_A, height[3], isfemale[2], isfemale[1], height[5]
internals = lp

Summary Statistics
    parameters mean std naive_se mcse ess rhat ess_per_sec 
        Symbol Float64 Float64 Float64 Float64 Float64 Float64 Float64 

   isfemale[1] 0.5061 0.5000 0.0050 0.0048 10062.6371 1.0003 15269.5556
   isfemale[2] 0.5108 0.4999 0.0050 0.0047 10069.3941 0.9999 15279.8090
   isfemale[3] 0.5116 0.4999 0.0050 0.0049 10003.7008 0.9999 15180.1226
   isfemale[4] 0.5114 0.4999 0.0050 0.0049 9756.3594 0.9999 14804.7942
   isfemale[5] 0.4988 0.5000 0.0050 0.0051 9872.9610 0.9999 14981.7314
   isfemale[6] 0.5141 0.4998 0.0050 0.0061 9637.4286 0.9999 14624.3226
   isfemale[7] 0.5134 0.4998 0.0050 0.0046 10231.4234 0.9999 15525.6804
   isfemale[8] 0.5025 0.5000 0.0050 0.0047 9834.8171 0.9999 14923.8499
   isfemale[9] 0.5144 0.4998 0.0050 0.0048 9842.9033 0.9999 14936.1203
  isfemale[10] 0.5200 0.4996 0.0050 0.0052 9767.5500 0.9999 14821.7754
     height[1] 170.4011 10.5277 0.1053 0.1129 9857.0370 1.0000 14957.5675
     height[2] 170.2644 10.7043 0.1070 0.0996 10133.5496 1.0000 15377.1618
     height[3] 170.4807 10.6361 0.1064 0.1049 9740.9282 1.0002 14781.3781
     height[4] 170.2575 10.6044 0.1060 0.1070 10072.0987 0.9999 15283.9130
     height[5] 170.5118 10.5506 0.1055 0.1080 9737.4231 0.9999 14776.0594
     height[6] 170.1812 10.6326 0.1063 0.1153 9745.2381 0.9999 14787.9181
     height[7] 170.3461 10.6142 0.1061 0.1094 9893.0560 1.0001 15012.2247
     height[8] 170.4455 10.5213 0.1052 0.0958 10029.7619 0.9999 15219.6690
     height[9] 170.4502 10.6084 0.1061 0.0981 9863.6789 0.9999 14967.6462
    height[10] 170.1680 10.5865 0.1059 0.1128 9526.2742 0.9999 14455.6513
      height_A 181.3947 4.6701 0.0467 0.0453 9569.4982 0.9999 14521.2416
    isfemale_A 0.0908 0.2873 0.0029 0.0028 9956.6954 0.9999 15108.7942

Quantiles
    parameters 2.5% 25.0% 50.0% 75.0% 97.5% 
        Symbol Float64 Float64 Float64 Float64 Float64 

   isfemale[1] 0.0000 0.0000 1.0000 1.0000 1.0000
   isfemale[2] 0.0000 0.0000 1.0000 1.0000 1.0000
   isfemale[3] 0.0000 0.0000 1.0000 1.0000 1.0000
   isfemale[4] 0.0000 0.0000 1.0000 1.0000 1.0000
   isfemale[5] 0.0000 0.0000 0.0000 1.0000 1.0000
   isfemale[6] 0.0000 0.0000 1.0000 1.0000 1.0000
   isfemale[7] 0.0000 0.0000 1.0000 1.0000 1.0000
   isfemale[8] 0.0000 0.0000 1.0000 1.0000 1.0000
   isfemale[9] 0.0000 0.0000 1.0000 1.0000 1.0000
  isfemale[10] 0.0000 0.0000 1.0000 1.0000 1.0000
     height[1] 151.0167 162.6419 170.2986 178.1220 190.5900
     height[2] 150.7162 162.2111 169.9460 178.2009 190.6453
     height[3] 150.9973 162.6603 170.1311 178.2281 190.7633
     height[4] 150.8729 162.4045 170.0385 178.1108 190.3861
     height[5] 151.2441 162.6010 170.2618 178.3460 190.7374
     height[6] 150.7536 162.2724 169.9650 177.9726 190.4018
     height[7] 150.4891 162.5470 170.2617 178.2245 190.2992
     height[8] 150.8853 162.7284 170.3147 178.1920 190.5763
     height[9] 151.0248 162.5023 170.1412 178.3231 190.8849
    height[10] 150.5369 162.2818 169.9072 177.9174 190.5773
      height_A 172.1184 178.2957 181.4849 184.5212 190.5233
    isfemale_A 0.0000 0.0000 0.0000 0.0000 1.0000

```

---

<div class="post-metadata">

**Author:** ![Hakan\_Kjellerstrand](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hakan_kjellerstrand/32/8447_2.png) [@Hakan\_Kjellerstrand](https://discourse.julialang.org/u/Hakan_Kjellerstrand)\
**Post date:** [July 16, 2021, 9:05pm UTC](https://discourse.julialang.org/t/help-with-turing-jl-model-using-categorical-normal-and-ordered-constraint/64765/3 "2021-07-16T21:05:06Z")

</div>

Thanks, @sethaxen .

I like your twist to skip the explicit ordering constraint and instead query about the order of the heights: `id_A = sortperm(height)[n-1]`

First I was a little confused that the listing of `height` has all about 170cm, but then understood what you do. And since I actually want the proper listing of all the (ordered) heights, I added `height_ordered`:

```julia
    height_ordered = Vector{T}(undef, n)
    perm = sortperm(height)
    for i in 1:n
        height_ordered[i] ~ Dirac(height[perm[n-i+1]])
    end

```

And now get this nice result:

```julia
   height_ordered[1] 186.5568 
   height_ordered[2] 181.2781 
   height_ordered[3] 177.6543 
   height_ordered[4] 174.4240 
   height_ordered[5] 171.4671 
   height_ordered[6] 168.5977 
   height_ordered[7] 165.7553 
   height_ordered[8] 162.7587 
   height_ordered[9] 159.3365 
  height_ordered[10] 154.5416

```

Again, thanks!