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
: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
Turing.@addlogprob! for observing the order, but there’s no difference between these two.
I guess that the real problem is the mix of
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?