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

I’m trying to model AllenDowney’s three problems “How tall is A?” (How tall is A? – Probably Overthinking It ):

``````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 ). Here’s the model of the second problem:

``````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 > height)
end

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

And it give the expected heights of:

``````   ...
height   176.4384    9.0135     0.2850    0.1777   824.6929    1.0007      179.5152
height   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):

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

``````   height   176.0685
height   184.1998
height   165.4559
height   161.9992
height   182.4204
height   175.6319
height   182.7167
height   161.5589
height   183.9556
height   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?

Here’s a model that does what you want:

``````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:

``````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, height, height, isfemale, isfemale, isfemale, isfemale, isfemale, height_A, isfemale, height, height, isfemale, isfemale, height, height, height, isfemale_A, height, isfemale, isfemale, height
internals         = lp

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

isfemale     0.5061    0.5000     0.0050    0.0048   10062.6371    1.0003    15269.5556
isfemale     0.5108    0.4999     0.0050    0.0047   10069.3941    0.9999    15279.8090
isfemale     0.5116    0.4999     0.0050    0.0049   10003.7008    0.9999    15180.1226
isfemale     0.5114    0.4999     0.0050    0.0049    9756.3594    0.9999    14804.7942
isfemale     0.4988    0.5000     0.0050    0.0051    9872.9610    0.9999    14981.7314
isfemale     0.5141    0.4998     0.0050    0.0061    9637.4286    0.9999    14624.3226
isfemale     0.5134    0.4998     0.0050    0.0046   10231.4234    0.9999    15525.6804
isfemale     0.5025    0.5000     0.0050    0.0047    9834.8171    0.9999    14923.8499
isfemale     0.5144    0.4998     0.0050    0.0048    9842.9033    0.9999    14936.1203
isfemale     0.5200    0.4996     0.0050    0.0052    9767.5500    0.9999    14821.7754
height   170.4011   10.5277     0.1053    0.1129    9857.0370    1.0000    14957.5675
height   170.2644   10.7043     0.1070    0.0996   10133.5496    1.0000    15377.1618
height   170.4807   10.6361     0.1064    0.1049    9740.9282    1.0002    14781.3781
height   170.2575   10.6044     0.1060    0.1070   10072.0987    0.9999    15283.9130
height   170.5118   10.5506     0.1055    0.1080    9737.4231    0.9999    14776.0594
height   170.1812   10.6326     0.1063    0.1153    9745.2381    0.9999    14787.9181
height   170.3461   10.6142     0.1061    0.1094    9893.0560    1.0001    15012.2247
height   170.4455   10.5213     0.1052    0.0958   10029.7619    0.9999    15219.6690
height   170.4502   10.6084     0.1061    0.0981    9863.6789    0.9999    14967.6462
height   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     0.0000     0.0000     1.0000     1.0000     1.0000
isfemale     0.0000     0.0000     1.0000     1.0000     1.0000
isfemale     0.0000     0.0000     1.0000     1.0000     1.0000
isfemale     0.0000     0.0000     1.0000     1.0000     1.0000
isfemale     0.0000     0.0000     0.0000     1.0000     1.0000
isfemale     0.0000     0.0000     1.0000     1.0000     1.0000
isfemale     0.0000     0.0000     1.0000     1.0000     1.0000
isfemale     0.0000     0.0000     1.0000     1.0000     1.0000
isfemale     0.0000     0.0000     1.0000     1.0000     1.0000
isfemale     0.0000     0.0000     1.0000     1.0000     1.0000
height   151.0167   162.6419   170.2986   178.1220   190.5900
height   150.7162   162.2111   169.9460   178.2009   190.6453
height   150.9973   162.6603   170.1311   178.2281   190.7633
height   150.8729   162.4045   170.0385   178.1108   190.3861
height   151.2441   162.6010   170.2618   178.3460   190.7374
height   150.7536   162.2724   169.9650   177.9726   190.4018
height   150.4891   162.5470   170.2617   178.2245   190.2992
height   150.8853   162.7284   170.3147   178.1920   190.5763
height   151.0248   162.5023   170.1412   178.3231   190.8849
height   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
``````
2 Likes

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`:

``````    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:

``````   height_ordered   186.5568
height_ordered   181.2781
height_ordered   177.6543
height_ordered   174.4240
height_ordered   171.4671
height_ordered   168.5977
height_ordered   165.7553
height_ordered   162.7587
height_ordered   159.3365
height_ordered   154.5416
``````

Again, thanks!