Unable to fit example Turing model

The following code, copied from the Turing website throws an error

using Turing
using StatsPlots

# Define a simple Normal model with unknown mean and variance.
@model gdemo(x, y) = begin
  s ~ InverseGamma(2, 3)
  m ~ Normal(0, sqrt(s))
  x ~ Normal(m, sqrt(s))
  y ~ Normal(m, sqrt(s))

#  Run sampler, collect results
chn = sample(gdemo(1.5, 2), HMC(0.1, 5), 1000)
ERROR: MethodError: no method matching AxisArrays.AxisArray(::Array{Float64,3}; iter
=1:1:1000, var=["m", "s", "lp", "n_steps", "is_accept", "acceptance_rate", "log_dens
ity", "hamiltonian_energy", "hamiltonian_energy_error", "step_size", "nom_step_size"
], chain=1:1)

I am on Julia 1.4.1 and Turing 0.10.1.

I duplicated an issue on the Turing GitHub here. The solution was to update AxisArrays to 0.4.2.