Continuing this thread…

A similar issue occurs with this PyMC3 model (also from Think Bayes):

```
k10 = 20 - 3
k01 = 15 - 3
k11 = 3
num_seen = k01 + k10 + k11
with pm.Model() as model9:
p0 = pm.Beta('p0', alpha=1, beta=1)
p1 = pm.Beta('p1', alpha=1, beta=1)
N = pm.DiscreteUniform('N', num_seen, 350)
q0 = 1-p0
q1 = 1-p1
ps = [q0*q1, q0*p1, p0*q1, p0*p1]
k00 = N - num_seen
data = pm.math.stack((k00, k01, k10, k11))
y = pm.Multinomial('y', n=N, p=ps, observed=data)
trace9 = pm.sample(1000, **options)
```

PyMC3 chose Metropolis:

```
Multiprocess sampling (2 chains in 2 jobs)
CompoundStep
>NUTS: [p1, p0]
>Metropolis: [N]
100.00% [4000/4000 00:02<00:00 Sampling 2 chains, 0 divergences]
Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 3 seconds.
The acceptance probability does not match the target. It is 0.5430480274605854, but should be close to 0.8. Try to increase the number of tuning steps.
The rhat statistic is larger than 1.05 for some parameters. This indicates slight problems during sampling.
The estimated number of effective samples is smaller than 200 for some parameters.
```

Now, to Julia:

```
@model function model9(k01, k10, k11)
p0 ~ Beta(1, 1)
p1 ~ Beta(1, 1)
num_seen = k01 + k10 + k11
N ~ DiscreteUniform(num_seen, 350)
q0 = 1-p0
q1 = 1-p1
ps = [q0*q1, q0*p1, p0*q1, p0*p1]
y = [N - num_seen, k01, k10, k11]
y ~ Multinomial(N, ps)
end
trace9 = sample(model9(k01, k10, k11), PG(50), 500)
```

We get:

```
Sampling: 100%|█████████████████████████████████████████| Time: 0:00:02
[112]:
Chains MCMC chain (500×9×1 Array{Float64, 3}):
Log evidence = 0.0
Iterations = 1:1:500
Number of chains = 1
Samples per chain = 500
Wall duration = 2.74 seconds
Compute duration = 2.74 seconds
parameters = p0, p1, N, y[1], y[2], y[3], y[4]
internals = lp, logevidence
Summary Statistics
parameters mean std naive_se mcse ess rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
p0 0.5122 0.2965 0.0133 0.0125 449.8449 0.9981 ⋯
p1 0.5232 0.2838 0.0127 0.0105 449.0627 0.9988 ⋯
N 190.5040 90.7350 4.0578 4.4045 514.1959 1.0049 ⋯
y[1] 44.0380 48.9169 2.1876 2.1886 572.9458 1.0044 ⋯
y[2] 47.2080 48.1557 2.1536 1.8543 469.3461 0.9985 ⋯
y[3] 47.9620 54.4876 2.4368 2.3048 503.4060 0.9983 ⋯
y[4] 51.2960 54.6575 2.4444 2.4123 430.4700 1.0021 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
p0 0.0350 0.2440 0.5235 0.7708 0.9749
p1 0.0447 0.2911 0.5222 0.7576 0.9844
N 39.9500 108.0000 191.5000 267.0000 338.5250
y[1] 0.0000 9.0000 25.0000 66.2500 180.7250
y[2] 0.4750 10.0000 30.0000 68.0000 183.0500
y[3] 0.0000 8.7500 27.0000 68.0000 199.2000
y[4] 0.0000 11.0000 32.0000 72.0000 185.6750
```

and:

This seems to be the same problem, but this time `PG(50)`

doesn’t work…

Neither does MH (which for some reason returns a normal distribution:

```
trace9 = sample(model9(k01, k10, k11), MH(), 500)
```