Unexpected error filling AxisArray with MvNormal

[markmbaum]

I'm trying to fill an AxisArray with values drawn from two MvNormal distributions from Distributions.jl. A snippet that should work is

using Random: Xoshiro, rand!
using Distributions
using AxisArrays

Xs = MvNormal([2000, 5000], [100 50; 50 100])
Xf = MvNormal([12000, 15000], [250 200; 200 250])
r = Xoshiro(10)
N = 100

zones = AxisArray(
    Matrix{Float64}(undef, 4, N),
    sample=[:sᵢ, :sₘ, :fᵢ, :fₘ],
    index=1:N
)

for i ∈ 1:N
    rand!(r, Xs, view(zones,1:2,i))
    rand!(r, Xf, view(zones,3:4,i))
end

but this raises an unexpected linear algebra error:

ERROR: MethodError: no method matching lmul!(::LinearAlgebra.LowerTriangular{Float64, LinearAlgebra.Adjoint{Float64, Matrix{Float64}}}, ::AxisArray{Float64, 1, SubArray{Float64, 1, Matrix{Float64}, Tuple{UnitRange{Int64}, Int64}, true}, Tuple{Axis{:sample, Vector{Symbol}}}})
Closest candidates are:
  lmul!(::Union{LinearAlgebra.LowerTriangular, LinearAlgebra.UnitLowerTriangular}, ::LinearAlgebra.LowerTriangular) at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/triangular.jl:1482
  lmul!(::LinearAlgebra.UniformScaling, ::AbstractVecOrMat) at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/uniformscaling.jl:305
  lmul!(::Number, ::AbstractArray) at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/generic.jl:219
  ...
Stacktrace:
 [1] unwhiten!(r::AxisArray{Float64, 1, SubArray{Float64, 1, Matrix{Float64}, Tuple{UnitRange{Int64}, Int64}, true}, Tuple{Axis{:sample, Vector{Symbol}}}}, a::PDMats.PDMat{Float64, Matrix{Float64}}, x::AxisArray{Float64, 1, SubArray{Float64, 1, Matrix{Float64}, Tuple{UnitRange{Int64}, Int64}, true}, Tuple{Axis{:sample, Vector{Symbol}}}})
   @ PDMats ~/.julia/packages/PDMats/CbBv1/src/generics.jl:42
 [2] unwhiten!(a::PDMats.PDMat{Float64, Matrix{Float64}}, x::AxisArray{Float64, 1, SubArray{Float64, 1, Matrix{Float64}, Tuple{UnitRange{Int64}, Int64}, true}, Tuple{Axis{:sample, Vector{Symbol}}}})
   @ PDMats ~/.julia/packages/PDMats/CbBv1/src/generics.jl:33
 [3] _rand!(rng::Xoshiro, d::FullNormal, x::AxisArray{Float64, 1, SubArray{Float64, 1, Matrix{Float64}, Tuple{UnitRange{Int64}, Int64}, true}, Tuple{Axis{:sample, Vector{Symbol}}}})
   @ Distributions ~/.julia/packages/Distributions/YQQXX/src/multivariate/mvnormal.jl:287
 [4] rand!(rng::Xoshiro, s::FullNormal, x::AxisArray{Float64, 1, SubArray{Float64, 1, Matrix{Float64}, Tuple{UnitRange{Int64}, Int64}, true}, Tuple{Axis{:sample, Vector{Symbol}}}})
   @ Distributions ~/.julia/packages/Distributions/YQQXX/src/genericrand.jl:91
 [5] top-level scope
   @ ~/Documents/code/marvin/mixing-monte-carlo/test.jl:32

although it does work if I fill the underlying array (zones.data) instead of the axis array. My versions are

AxisArrays v0.4.6
Distributions v0.25.86