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[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:
...
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):
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[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?