Dispatch error on loglikelihood when sampling from arraydist of multivariate distributions

[PTWaade]

Minimal working example

using Turing

@model function m1(x) #will error

    n ~ Normal(0, 1)
    
    dists = [
        MvNormal([n,n], I),
        MvNormal([n,n], I),
        MvNormal([n,n], I),
        ]

    x ~ arraydist(dists)
end

@model function m2(x) #will work

    n ~ Normal(0, 1)
    
    dists = [
        MvNormal([n,n], I),
        MvNormal([n,n], I),
        MvNormal([n,n], I),
        ]

    for i in length(x)
        x[i] ~ dists[i]
    end
end


@model function m3() #will work

    n ~ Normal(0, 1)
    
    dists = [
        MvNormal([n,n], I),
        MvNormal([n,n], I),
        MvNormal([n,n], I),
        ]

    x ~ arraydist(dists)
end


model = m1([[1., 1.], [1., 1.], [1., 1.]])
sample(model, NUTS(), 1000)

model = m2([[1., 1.], [1., 1.], [1., 1.]])
sample(model, NUTS(), 1000)

model = m3()
sample(model, NUTS(), 1000)

Description

Hey Turing,

As always, a huge appreciation for the incredible software you are building!

I hope this is not a known issue or a duplicate!

I get this error:

MethodError: no method matching loglikelihood(::DistributionsAD.VectorOfMultivariate{Continuous, IsoNormal, Vector{…}}, ::Vector{Vector{…}})

when I sample from an arraydist of multivariate distributions (see model m1).

It works if I sample in a for loop (m2), and it also works when there is no data and it estimates x as parameters (m3). So I expect a fix shouldn't be too difficult, although I don't know the internals of Turing's calculation of loglikelihoods well enough to say for sure.

Let me know if I can do anything! The easy workaround for me is of course just to use the for loop for now.

Warm regards!

Julia version info

versioninfo()

Julia Version 1.11.5
Commit 760b2e5b739 (2025-04-14 06:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin24.0.0)
  CPU: 8 × Apple M1
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, apple-m1)
Threads: 2 default, 0 interactive, 1 GC (on 4 virtual cores)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 2

Manifest

]st --manifest

[47edcb42] ADTypes v1.14.0
  [621f4979] AbstractFFTs v1.5.0
  [80f14c24] AbstractMCMC v5.6.0
  [7a57a42e] AbstractPPL v0.11.0
  [1520ce14] AbstractTrees v0.4.5
  [7d9f7c33] Accessors v0.1.42
  [79e6a3ab] Adapt v4.3.0
  [0bf59076] AdvancedHMC v0.7.1
  [5b7e9947] AdvancedMH v0.8.7
⌃ [576499cb] AdvancedPS v0.6.1
⌅ [b5ca4192] AdvancedVI v0.2.11
  [66dad0bd] AliasTables v1.1.3
  [dce04be8] ArgCheck v2.5.0
  [4fba245c] ArrayInterface v7.18.1
  [a9b6321e] Atomix v1.1.1
  [13072b0f] AxisAlgorithms v1.1.0
  [39de3d68] AxisArrays v0.4.7
  [198e06fe] BangBang v0.4.4
  [9718e550] Baselet v0.1.1
  [76274a88] Bijectors v0.15.6
  [082447d4] ChainRules v1.72.3
  [d360d2e6] ChainRulesCore v1.25.1
  [0ca39b1e] Chairmarks v1.3.1
  [9e997f8a] ChangesOfVariables v0.1.10
  [861a8166] Combinatorics v1.0.3
  [38540f10] CommonSolve v0.2.4
  [bbf7d656] CommonSubexpressions v0.3.1
  [34da2185] Compat v4.16.0
  [a33af91c] CompositionsBase v0.1.2
  [88cd18e8] ConsoleProgressMonitor v0.1.2
  [187b0558] ConstructionBase v1.5.8
  [a8cc5b0e] Crayons v4.1.1
  [9a962f9c] DataAPI v1.16.0
  [864edb3b] DataStructures v0.18.22
  [e2d170a0] DataValueInterfaces v1.0.0
  [244e2a9f] DefineSingletons v0.1.2
  [8bb1440f] DelimitedFiles v1.9.1
  [b429d917] DensityInterface v0.4.0
  [163ba53b] DiffResults v1.1.0
  [b552c78f] DiffRules v1.15.1
  [a0c0ee7d] DifferentiationInterface v0.6.53
  [31c24e10] Distributions v0.25.119
  [ced4e74d] DistributionsAD v0.6.58
  [ffbed154] DocStringExtensions v0.9.4
⌃ [366bfd00] DynamicPPL v0.36.1
  [cad2338a] EllipticalSliceSampling v2.0.0
  [4e289a0a] EnumX v1.0.5
  [e2ba6199] ExprTools v0.1.10
  [55351af7] ExproniconLite v0.10.14
  [7a1cc6ca] FFTW v1.8.1
  [9aa1b823] FastClosures v0.3.2
  [1a297f60] FillArrays v1.13.0
  [6a86dc24] FiniteDiff v2.27.0
⌅ [f6369f11] ForwardDiff v0.10.38
  [069b7b12] FunctionWrappers v1.1.3
  [77dc65aa] FunctionWrappersWrappers v0.1.3
  [d9f16b24] Functors v0.5.2
  [46192b85] GPUArraysCore v0.2.0
  [076d061b] HashArrayMappedTries v0.2.0
  [34004b35] HypergeometricFunctions v0.3.28
  [22cec73e] InitialValues v0.3.1
  [a98d9a8b] Interpolations v0.15.1
  [8197267c] IntervalSets v0.7.11
  [3587e190] InverseFunctions v0.1.17
  [41ab1584] InvertedIndices v1.3.1
  [92d709cd] IrrationalConstants v0.2.4
  [c8e1da08] IterTools v1.10.0
  [82899510] IteratorInterfaceExtensions v1.0.0
  [692b3bcd] JLLWrappers v1.7.0
  [682c06a0] JSON v0.21.4
  [ae98c720] Jieko v0.2.1
  [63c18a36] KernelAbstractions v0.9.34
  [5ab0869b] KernelDensity v0.6.9
  [5be7bae1] LBFGSB v0.4.1
  [8ac3fa9e] LRUCache v1.6.2
  [b964fa9f] LaTeXStrings v1.4.0
  [1d6d02ad] LeftChildRightSiblingTrees v0.2.0
⌅ [6f1fad26] Libtask v0.8.8
  [d3d80556] LineSearches v7.3.0
  [6fdf6af0] LogDensityProblems v2.1.2
  [996a588d] LogDensityProblemsAD v1.13.0
  [2ab3a3ac] LogExpFunctions v0.3.29
  [e6f89c97] LoggingExtras v1.1.0
  [c7f686f2] MCMCChains v6.0.7
  [be115224] MCMCDiagnosticTools v0.3.14
  [e80e1ace] MLJModelInterface v1.11.1
  [1914dd2f] MacroTools v0.5.16
  [dbb5928d] MappedArrays v0.4.2
  [128add7d] MicroCollections v0.2.0
  [e1d29d7a] Missings v1.2.0
  [2e0e35c7] Moshi v0.3.5
  [d41bc354] NLSolversBase v7.9.1
  [872c559c] NNlib v0.9.30
  [77ba4419] NaNMath v1.1.3
  [86f7a689] NamedArrays v0.10.3
  [c020b1a1] NaturalSort v1.0.0
  [6fe1bfb0] OffsetArrays v1.17.0
  [429524aa] Optim v1.12.0
  [3bd65402] Optimisers v0.4.6
  [7f7a1694] Optimization v4.2.0
  [bca83a33] OptimizationBase v2.5.0
  [36348300] OptimizationOptimJL v0.4.3
  [bac558e1] OrderedCollections v1.8.0
  [90014a1f] PDMats v0.11.34
  [d96e819e] Parameters v0.12.3
  [69de0a69] Parsers v2.8.3
  [85a6dd25] PositiveFactorizations v0.2.4
⌅ [aea7be01] PrecompileTools v1.2.1
  [21216c6a] Preferences v1.4.3
  [08abe8d2] PrettyTables v2.4.0
  [33c8b6b6] ProgressLogging v0.1.4
  [92933f4c] ProgressMeter v1.10.4
  [43287f4e] PtrArrays v1.3.0
  [1fd47b50] QuadGK v2.11.2
  [74087812] Random123 v1.7.0
  [e6cf234a] RandomNumbers v1.6.0
  [b3c3ace0] RangeArrays v0.3.2
  [c84ed2f1] Ratios v0.4.5
  [c1ae055f] RealDot v0.1.0
  [3cdcf5f2] RecipesBase v1.3.4
  [731186ca] RecursiveArrayTools v3.33.0
  [189a3867] Reexport v1.2.2
  [ae029012] Requires v1.3.1
  [79098fc4] Rmath v0.8.0
  [f2b01f46] Roots v2.2.7
  [7e49a35a] RuntimeGeneratedFunctions v0.5.14
⌅ [26aad666] SSMProblems v0.1.1
  [0bca4576] SciMLBase v2.88.0
  [c0aeaf25] SciMLOperators v0.4.0
  [53ae85a6] SciMLStructures v1.7.0
  [30f210dd] ScientificTypesBase v3.0.0
  [7e506255] ScopedValues v1.3.0
  [efcf1570] Setfield v1.1.2
  [a2af1166] SortingAlgorithms v1.2.1
  [9f842d2f] SparseConnectivityTracer v0.6.18
  [dc90abb0] SparseInverseSubset v0.1.2
  [0a514795] SparseMatrixColorings v0.4.19
  [276daf66] SpecialFunctions v2.5.1
  [171d559e] SplittablesBase v0.1.15
  [90137ffa] StaticArrays v1.9.13
  [1e83bf80] StaticArraysCore v1.4.3
  [64bff920] StatisticalTraits v3.4.0
  [10745b16] Statistics v1.11.1
  [82ae8749] StatsAPI v1.7.0
  [2913bbd2] StatsBase v0.34.5
  [4c63d2b9] StatsFuns v1.5.0
  [892a3eda] StringManipulation v0.4.1
  [09ab397b] StructArrays v0.7.1
  [2efcf032] SymbolicIndexingInterface v0.3.40
  [3783bdb8] TableTraits v1.0.1
  [bd369af6] Tables v1.12.0
  [5d786b92] TerminalLoggers v0.1.7
  [9f7883ad] Tracker v0.2.38
  [28d57a85] Transducers v0.4.84
  [fce5fe82] Turing v0.38.0
  [3a884ed6] UnPack v1.0.2
  [013be700] UnsafeAtomics v0.3.0
  [efce3f68] WoodburyMatrices v1.0.0
  [700de1a5] ZygoteRules v0.2.7
  [f5851436] FFTW_jll v3.3.11+0
  [1d5cc7b8] IntelOpenMP_jll v2025.0.4+0
  [81d17ec3] L_BFGS_B_jll v3.0.1+0
  [856f044c] MKL_jll v2025.0.1+1
  [efe28fd5] OpenSpecFun_jll v0.5.6+0
  [f50d1b31] Rmath_jll v0.5.1+0
  [1317d2d5] oneTBB_jll v2022.0.0+0
  [0dad84c5] ArgTools v1.1.2
  [56f22d72] Artifacts v1.11.0
  [2a0f44e3] Base64 v1.11.0
  [ade2ca70] Dates v1.11.0
  [8ba89e20] Distributed v1.11.0
  [f43a241f] Downloads v1.6.0
  [7b1f6079] FileWatching v1.11.0
  [9fa8497b] Future v1.11.0
  [b77e0a4c] InteractiveUtils v1.11.0
  [4af54fe1] LazyArtifacts v1.11.0
  [b27032c2] LibCURL v0.6.4
  [76f85450] LibGit2 v1.11.0
  [8f399da3] Libdl v1.11.0
  [37e2e46d] LinearAlgebra v1.11.0
  [56ddb016] Logging v1.11.0
  [d6f4376e] Markdown v1.11.0
  [a63ad114] Mmap v1.11.0
  [ca575930] NetworkOptions v1.2.0
  [44cfe95a] Pkg v1.11.0
  [de0858da] Printf v1.11.0
  [3fa0cd96] REPL v1.11.0
  [9a3f8284] Random v1.11.0
  [ea8e919c] SHA v0.7.0
  [9e88b42a] Serialization v1.11.0
  [1a1011a3] SharedArrays v1.11.0
  [6462fe0b] Sockets v1.11.0
  [2f01184e] SparseArrays v1.11.0
  [f489334b] StyledStrings v1.11.0
  [4607b0f0] SuiteSparse
  [fa267f1f] TOML v1.0.3
  [a4e569a6] Tar v1.10.0
  [8dfed614] Test v1.11.0
  [cf7118a7] UUIDs v1.11.0
  [4ec0a83e] Unicode v1.11.0
  [e66e0078] CompilerSupportLibraries_jll v1.1.1+0
  [deac9b47] LibCURL_jll v8.6.0+0
  [e37daf67] LibGit2_jll v1.7.2+0
  [29816b5a] LibSSH2_jll v1.11.0+1
  [c8ffd9c3] MbedTLS_jll v2.28.6+0
  [14a3606d] MozillaCACerts_jll v2023.12.12
  [4536629a] OpenBLAS_jll v0.3.27+1
  [05823500] OpenLibm_jll v0.8.5+0
  [bea87d4a] SuiteSparse_jll v7.7.0+0
  [83775a58] Zlib_jll v1.2.13+1
  [8e850b90] libblastrampoline_jll v5.11.0+0
  [8e850ede] nghttp2_jll v1.59.0+0
  [3f19e933] p7zip_jll v17.4.0+2

[penelopeysm]

This comes from the use of loglikelihood in observe here:

https://github.com/TuringLang/DynamicPPL.jl/blob/cdeb65752297009ef1f2e2eef95e3b8b963e5cff/src/context_implementations.jl#L290-L293

loglikelihood isn't implemented for filldist and arraydist, so any observed variable with these distributions will fail.

On the other hand, assumed variables (i.e. latent parameters) call logpdf instead of loglikelihood, which is why m3() doesn't error.

https://github.com/TuringLang/DynamicPPL.jl/blob/cdeb65752297009ef1f2e2eef95e3b8b963e5cff/src/context_implementations.jl#L284-L285

The cheapest fix would be to patch DistributionsAD to include loglikelihood methods for those. We already have logpdf calculations in there so I suspect it would just be a few lines of changes

• https://github.com/TuringLang/DistributionsAD.jl/blob/main/src/arraydist.jl
• https://github.com/TuringLang/DistributionsAD.jl/blob/main/src/filldist.jl

But more generally, I don't fully understand why we use loglikelihood for observe, but logpdf for assume, it can lead to inconsistent behaviour like this. @mhauru and I had a little gripe about it before, but I'm not sure if we actually had a reason why we couldn't use logpdf for both?

[PTWaade]

I realised (feeling a little foolish) that it works fine when just entering x as a matrix instead of a vector of vectors. This works:

model = m1([[1. 1. 1. ; 1. 1. 1.]])
sample(model, NUTS(), 1000)

Maybe one could automate the transformation there, or have a nice error message that saves people like me from making this mistake (using a vector of vectors feels very intuitive).

I closed the issue (in embarrassment), but saw that you responded with a somewhat more general point - feel free to reopen!

And thanks for the fast response :)

[PTWaade]

(The vector of vectors solution would also be more flexible, since it wouldn't require all the vectors to be the same length, so perhaps this would be nice to support after all)

[PTWaade]

Sorry for the many messages: I just realised it of course should have been:

model = m1([1. 1. 1. ; 1. 1. 1.])
sample(model, NUTS(), 1000)

Without wrapping the matrix in a vector. Both works, however - and it surprises me quite a lot that the matrix-wrapped-in-vector option worked

[mhauru]

Glad you figured out how to get it to work. I don't think we'll add support for vectors of vectors. I'm not sure arraydist supports different dimensions for the different column distributions either (or is it row distributions? I always mix this up.). A better error message could be nice tough, maybe in Distributions.jl, saying "we expected a matrix but you gave us a vector of vectors".

I, too, was surprised that the matrix wrapped in a vector worked. I suspect the answer is automatic broadcasting for loglikelihood.

But more generally, I don't fully understand why we use loglikelihood for observe, but logpdf for assume, it can lead to inconsistent behaviour like this. @mhauru and I had a little gripe about it before, but I'm not sure if we actually had a reason why we couldn't use logpdf for both?

I tried switching to logpdf everywhere a couple of months ago, the problem I ran into was related to the above implicit broadcasting. We allow this:

julia> @model function f(vector, matrix)
           mean ~ Normal()
           vector ~ Normal(mean)
           matrix ~ MvNormal(fill(mean, 2), I)
           return nothing
       end

julia> f([1.0, 2.0], [1.0 2.0 3.0; 4.0 5.0 6.0])

logpdf and loglikelihood handle this sort of implicit broadcasting differently:

julia> logpdf(Normal(), [0.1, 0.2])
┌ Warning: `logpdf(d::UnivariateDistribution, X::AbstractArray{<:Real})` is deprecated, use `map(Base.Fix1(logpdf, d), X)` instead.
│   caller = top-level scope at REPL[23]:1
└ @ Core REPL[23]:1
2-element Vector{Float64}:
 -0.9239385332046728
 -0.9389385332046728

julia> loglikelihood(Normal(), [0.1, 0.2])
-1.8628770664093457

and we need the latter behaviour. This only comes up with observe because with assume the variable can always be assumed to be of the dimension of the distribution, since it's never provided by the user.

[penelopeysm]

Oh, yes, but I was thinking we should disallow it.

[mhauru]

I think we would be better off without it (since .~ exists), I guess the question is whether it's worth the breakage. Opened a discussion about it here.