Unexpected broadcasting behavior of pdf function

[OkonSamuel]

@bkamins recently pointed out attention to an unexpected behavior that occurs when broadcasting with the pdf function. This issue tracks it.

julia> using CategoricalDistributions

julia> import CategoricalDistributions: classes

julia> u = UnivariateFinite(['x', 'z'], rand(2, 3, 2), pool=missing, ordered=true)
2×3 UnivariateFiniteArray{OrderedFactor{2}, Char, UInt8, Float64, 2}:
 UnivariateFinite{OrderedFactor{2}}(x=>0.215, z=>0.742)  …  UnivariateFinite{OrderedFactor{2}}(x=>0.903, z=>0.735)
 UnivariateFinite{OrderedFactor{2}}(x=>0.575, z=>0.907)     UnivariateFinite{OrderedFactor{2}}(x=>0.911, z=>0.904)

julia> l = u[1:2, 1:1]
2×1 UnivariateFiniteArray{OrderedFactor{2}, Char, UInt8, Float64, 2}:
 UnivariateFinite{OrderedFactor{2}}(x=>0.215, z=>0.742)
 UnivariateFinite{OrderedFactor{2}}(x=>0.575, z=>0.907)

julia> c = permutedims(classes(l))
1×2 CategoricalArrays.CategoricalArray{Char,2,UInt8}:
 'x'  'z'

julia> pdf.(l, c) # Expected a (2,2) matrix
2×1 Matrix{Float64}:
 0.21520308715626935
 0.9066550037643591

As shown in the above code block, the output of pdf(l,c) was expected 2 x 2 matrix, but instead gave a 2 x 1 matrix.
The relevant section of our codebase is https://github.com/JuliaAI/CategoricalDistributions.jl/blob/dev/src/arrays.jl#L184-L204.

[OkonSamuel]

On further investigation into the issue. I discovered that the existing broadcasting functionality defined here only works as expected when length.(axes(u)) == length.(axes(v)). For example, the code block above gave unexpected output because length.(axes(l)) == length.(axes(c)) is false

julia> l
2×1 UnivariateFiniteArray{OrderedFactor{2}, Char, UInt8, Float64, 2}:
 UnivariateFinite{OrderedFactor{2}}(x=>0.215, z=>0.742)
 UnivariateFinite{OrderedFactor{2}}(x=>0.575, z=>0.907)

julia> c
1×2 CategoricalArrays.CategoricalArray{Char,2,UInt8}:
 'x'  'z'

julia> length.(axes(l)) == length.(axes(c))
false