Accept Number types in distribution constructors.

[rafaqz]

There are many Number typed objects that would work just fine in a distribution but currently don't. #1413 has examples where units would already work in a lot of cases by just using Number.

ModelParameters.jl is a similar case. It has Number typed parameters Param so that e.g. Unitful.jl units can be put in them.

But then they don't work as input to Distributions.jl objects when users try to do that:
rafaqz/ModelParameters.jl#51

Embedding parametrised distributions in large nested models could be extremely easy with this method.

The alternative is we add a parallel RealParam object and add new types and methods to AbstractNumbers.jl to work on Real. This is a lot of work and code duplication compared to simply swapping Real to Number here.

[devmotion]

I thought #1413 showed that just replacing Real with Number would lead to many functions returning incorrect units?

[rafaqz]

That's not the impression I got - there are some errors, but a lot of things actually just work? Like these examples #1413 (comment)

I guess the question is do you want to guarantee that there are never any wrong outputs. And accepting Number means you cant. I get why that's a problem.

But it leaves this impasse where it is impossible to generically wrap numbers in another object and let them be useful in the ecosystem.

We really need a trait-based number system in base so you could call isreal, iscomplex on a number rather than being stuck in the hierarchy.

Then these could be checked in the constructor to guarantee correct results, rather than fixing the type to Real.

[mschauer]

Relatedly, for example, https://github.com/gaurav-arya/StochasticAD.jl/blob/576d4e9026e2b4c225b2aab5910f4b6148a66666/src/discrete_randomness.jl#L110 is to circumvent that the Binomials don't allow numbers @gaurav-arya

[devmotion]

Regarding the example, I think it's annoying if constructors "lie" about the type of the instance they return.

I think it would be nice to generalize the definitions. I'm just not a fan of simple search-replace approaches, I think it deserves to take one distribution, put together a PR for this specific distribution only, and discuss it intensively to figure out all problems and fix them. That would make the overall design and potential pitfalls much clearer. Afterwards one could generalize it to other distributions.

Another caveat, as you mention, is that generally we really don't want arbitrary but only real numbers as distribution parameters and arguments to e.g. pdf and cdf (one exception would be cf and mgf: #1504). To prevent issues with user-provided inputs, we could check the arguments and parameters with isreal.

[rafaqz]

we could check the arguments and parameters with isreal

This sounds like a solution to me.

Can we switch the type to Number but have an an isreal check in the constructors? rather than limiting the type to Real ?

ModelParameters.jl Param has no numerical properties different from the wrapped number, except that we know it is the parameter of a model, and we can find in any constructed object.

So there are no meaningful PRs to make except this change to Number instead of Real, and checking with isreal.

[rafaqz]

Also I may not have communicated clearly enough: The main ModelParameters.jl use case here is not with unitful units. Its just with Real.

But to also allow units in other contexts, out Param wrappers have to be Number. So we are stuck between the competing needs of two packages.

[devmotion]

On a second thought, I assume one shouldn't expect too much from loosening types in Distributions - for many distributions, evaluations are based on functions in StatsFuns or SpecialFunctions which often only support Real. Many methods are even highly optimized/only exist for Float64 and possibly Float16/Float32. Even Dual numbers are a problem sometimes and require custom forward mode rules. So it's not completely obvious for which distributions changes in Distributions would be sufficient.

In my experience wrapper types lead to all kind of problems and are usually less convenient to work with than the unwrapped values, even if one tries to widen function signatures. Could you automatically unwrap the parameters when constructing distributions?

[mschauer]

Unwrapping the parameter might introduce the need to wrap the distribution object which then causes the annoying constructor dishonesty (i. E instead of a distribution with dual parameters you get a wrapper type with a distribution and it’s dual distribution counterpart)

[gaurav-arya]

A quick clarification: what we need in StochasticAD is to relax Int to Integer, rather than reals to numbers, which is a bit different and I would guess easier to do.

[rafaqz]

ModelParameters.jl Param wrappers don't lead to any of these problems, and the act of wrapping is the whole and only point of the package.

They are usually stripped out of a model before use. And if not, they wont actually flow through other code after being stored in the distribution - any operation on them will be on the parent number.

They just provide a way to build arbitrary nested model structures but update the parameters in an automated way.

For example, by putting a Param in this distribution:

somestructfield = Exponential(Param(1.0))

it can be found and rebuilt as

somestructfield = Exponential(1.0)

or when e.g. an optimiser wants to update it, to:

somestructfield = Exponential(0.9)

But currently this is not possible, because they are Numbers. Loosening the type will absolutely solve this problem.

[devmotion]

I'll have a look at the package before commenting any further on what ModelParameters does 🙂 I've only worked with https://github.com/invenia/ParameterHandling.jl a few times which seems to handle parameters and their optimization by defining a function that constructs the model from the given parameters.

[rafaqz]

Ok thanks, thats a good idea. A good place to start is this issue from a user of both of these packages:
rafaqz/ModelParameters.jl#51

Op essentially sums up why ModelParameters.jl is useful to him and why it should work with Distributions.jl. I didn't quite understand him until I tried to do the same myself. But it would be very elegant to put a Param in a distribution in a model definition.

You will also notice that ParameterHandling.jl jumps through some hoops do things what would be very simple using a Param as I outline above:

deferred has some type-stability issues when used in conjunction with abstract types. For example, flatten(deferred(Normal, 5.0, 4.0)) won't infer properly. A simple work around is to write a function normal(args...) = Normal(args...) and work with deferred(normal, 5.0, 4.0) instead.

[devmotion]

I read through the issue and the README of the package, and the use case is a bit clearer to me now. However, to me it seems that generally ModelParameters.jl faces the same limitations as e.g. ForwardDiff by design: you can only include parameterized components in the (sub)models if the corresponding structs are written generically enough and support Params. That's a quite strong assumption in general, even more so since Params is not a subtype of Real, in contrast to ForwardDiff.Dual.

If you just define a function that constructs your model from a vector/named tuple/... of parameters (manually or e.g. with ParameterHandling), you don't run into such issues. As far as I can tell, the mentioned problem with deferred is orthogonal to this approach: one can define model building functions without using deferred (im- or explicitly) and this type inference issue occurs in other use cases as well (see e.g. TuringLang/Turing.jl#1934 (comment)).

[rafaqz]

Are you saying you don't want to relax the type to Number and check the input for isreal?

(You seem to be missing that the type is Number to also allow Unitful.jl and many other similar packages that use Number - this is not a random decision made in a vacuum just for my package - Unitful.jl support is much more important for physical and process based modelling than Distributions.jl support. People just want both)

[devmotion]

I'm saying that I think generally you can't expect it to "just work" but this approach requires you to change every package you want to use if it is not written generically enough.

Regarding Distributions, I think we could try to change Real to Number for one example distribution (e.g., Normal) and see how many changes are needed and what the resulting implementation would look like, so we can decide if we want to go this route. But even in this simple case we already know that it won't be sufficient to change the code in Distributions: pdf etc. won't work for parameters of type Number because that is not supported by StatsFuns (https://github.com/JuliaStats/StatsFuns.jl/blob/master/src/distrs/norm.jl). In this case we could update StatsFuns as well but for many other distributions we rely on algorithms in SpecialFunctions for which these updates are not possible or at least more challenging.

Edit: It's even worse, I missed that even for Normal we use SpecialFunctions: the cdf uses erfc which is only implemented for FloatXX and BigFloat (and inputs x where float(x) is one of these types).

[rafaqz]

Ok, I'll play with doing that with Normal and make a PR.

The other issues may affect other Number types, but ModelParameters.jl Param can be stripped out of the model automatically before use.

The problem is just being able to construct distribution objects with Param wrappers in them.

[rafaqz]

Of course the other option is that we write a parallel RealParam wrapper in ModelParameters.jl (it has to separately implement the whole Real number interface, so not as simple as it sounds)

It just feels like a flaw in the design of Julias number system that we have to do that kind of duplication. There should have been more reliance on traits.

[gaurav-arya]

@rafaqz regarding the prospect of writing a RealParam -- I have also been working with overloading generic number types and thinking about this sort of thing, although I haven't actually done it yet.

Obviously, we're up against a Julia type system limitation, which is that we can't do MyType{T} <: T. But I somewhat agree with the sentiment that given this limitation, the "onus" is somewhat on the developers of these number types to hack around the system to get this to work, and that's the only general way to do it in current Julia with subtype-based dispatches abound. (Although, the Unitful argument is a separate thing, so not commenting on the calculus for making Distributions.jl specifically more accepting.)

I do wonder how ergonomic this hacking could be for a specific target, e.g. Real though, using abstract typing of your generic number types to avoid code dup as well as liberal use of things like @forward. I would guess what's done in Symbolics.jl might be the state of the art here, but I'm also still learning. And also, if it really holds that number types could pass through a lot of Distributions.jl, then you could initially be lazy about supporting a lot of the real-number interface for RealParam and things would work fine without at least having to fight the dispatch system.

[rafaqz]

Yeah, its not the easiest situation.

AbstractNumbers.jl was the easiest fix for me, although I'm the defacto maintainer these days:
https://github.com/SimonDanisch/AbstractNumbers.jl

I use it to make Param work as a basic number. Though it can be stripped out so easily its best if it still basically works and handles any checks that are done in constructors.

My idea was to just add a new AbstractNumbers.AbstractReal <: Real type to AbstractNumbers.jl and switch dispatch of all the methods it extends to Union{AbstractReal,AbstractNumber}

Then have RealParam <: AbstractNumbers.AbstractReal and NumberParam <: AbstractNumbers.AbstractNumber and use a par(x) function in ModelParameters.jl to choose them based on the type of x.

I'm not sure what Symbolics.jl does.

[DNF2]

The constructors seem very inconsistent to me:

julia> Uniform(0, 1)
Uniform{Float64}(a=0.0, b=1.0)

julia> Uniform{Float64}(0, 1)
Uniform{Float64}(a=0.0, b=1.0)

julia> Normal(0, 1)
Normal{Float64}(μ=0.0, σ=1.0)

julia> Normal{Float64}(0, 1)
ERROR: MethodError: no method matching Normal{Float64}(::Int64, ::Int64)

This is all over the place. Also, there are unnecessary constructors, like

struct Normal{T<:Real} <: ContinuousUnivariateDistribution
    μ::T
    σ::T
    Normal{T}(µ::T, σ::T) where {T<:Real} = new{T}(µ, σ)  # <= why not use default constructor?
end

There are many constructors that seem redundant, where it would be better to leave it to default constructors.

I also wonder why the @check_args is happening in outer constructors, instead of inner. So you have

julia> Normal(0.0, -1.0)
ERROR: DomainError with -1.0:

But this is fine:

julia> Normal{Float64}(0.0, -1.0)
Normal{Float64}(μ=0.0, σ=-1.0)

This one is strange to me:

julia> Uniform(-1, 1)
Uniform{Float64}(a=-1.0, b=1.0)

julia> Uniform{Int}(-1, 1)
Uniform{Int64}(a=-1, b=1)

If Integer parameters is ok, why auto-convert to float?

[devmotion]

The constructors seem very inconsistent to me

Unfortunately, many things in Distributions are not completely consistent. My impression is that this is mainly due to how the package evolved: it's an old package and older Julia versions required different syntax and implementations. So some things were just not possible to do in the initial stages and some design choices only turned out to be suboptimal later.

I also wonder why the @check_args is happening in outer constructors, instead of inner.

These checks cause overhead (noticeable eg in AD which required some special fixes) and hence it is nice to be able to skip them if not needed (e.g. in convert or constructors with default parameters).

If Integer parameters is ok, why auto-convert to float?

This is mainly for historical reasons. It took quite some effort to make pdf etc type stable and without this conversion it was quite easy to run into these type stability issues before.

[DNF2]

My impression is that this is mainly due to how the package evolved: it's an old package and older Julia versions required different syntax and implementations.

With respect to the current issue, if there is any attempt to widen the input types to Number, then harmonizing the constructors could also be considered. In some cases it would be as easy as just removing the code and relying on defaults.

These checks cause overhead

In that case check_args could be explicitly called with false? Some extra verbosity could be ok for performance optimization, espcially since it is not safe.

But the worst inconsistency is that some have constructors

A{T}(µ::T, σ::T) where {T<:Real} = new{T}(µ, σ)

and some have

B{T}(µ, σ) where {T<:Real} = new{T}(µ, σ)

In practice, they could both simply be deleted.

[devmotion]

In that case check_args could be explicitly called with false?

Just passing false was not sufficient to fix performance issues for AD. Hence I became convinced that it's desirable to not perform checks in inner constructors - and that it's also reasonable. Users are supposed to use outer constructors which promote arguments etc. - whereas if you call an inner constructor, you specify the desired type explicitly and hence indicate that you really want an instance of some specific type.

So IMO we want inner constructors without checks + outer constructors with checks (if needed).

[DNF2]

If it's a deliberate design decision, I won't quibble. I could imagine a different interface, though. For example Normal{T}(a, b, NoArgCheck()) etc.