Float32

.gitignore

@@ -11,3 +11,4 @@ docs/Manifest.toml
 /jmd/*.md
 /jmd/figures/
 
+test/Manifest.toml

docs/make.jl

@@ -20,13 +20,14 @@ makedocs(;
     ),
     pages=[
         "Home" => "index.md",
-        "Dependencies" => "set_optimize.md",
+        "Parameter type" => "partype.md",
         "Distributions" => [
             "LogNormal" => "lognormal.md",
             "LogitNormal" => "logitnormal.md",
             "Weibull" => "weibull.md",
             "Gamma" => "gamma.md",
         ],
+        "Dependencies" => "set_optimize.md",
         "API" => "api.md",
         #"Details" => "z_autodocs.md",
     ],

docs/src/index.md

@@ -20,7 +20,7 @@ d = fit(LogNormal, Moments(3.0,4.0))
 ```
 - mean and upper quantile point
 ```jldoctest; output = false
-d = fit(LogNormal, 3, @qp_uu(8))
+d = fit(LogNormal, 3.0, @qp_uu(8))
 (mean(d), quantile(d, 0.975)) .≈ (3.0, 8.0)
 # output
 (true, true)
@@ -30,7 +30,7 @@ d = fit(LogNormal, 3, @qp_uu(8))
 DocTestSetup = :(using Statistics,DistributionFits,Optim)
 ```
 ```jldoctest; output = false
-d = fit(LogNormal, 3, @qp_uu(8), Val(:mode))
+d = fit(LogNormal, 3.0, @qp_uu(8), Val(:mode))
 (mode(d), quantile(d, 0.975)) .≈ (3.0, 8.0)
 # output
 (true, true)
@@ -73,10 +73,6 @@ d = fit(LogNormal, moments(dn))
 fit(::Type{D}, ::AbstractMoments) where {D<:Distribution}
 ```
 
-```@docs
-moments(d::Distribution, ::Val{N} = Val(2)) where N 
-```
-
 ## Fit to several quantile points
 
 ```@docs

docs/src/lognormal.md

@@ -37,7 +37,7 @@ true
 ```
 
 ```@docs
-fit(::Type{LogNormal}, ::Any, ::AbstractΣstar) 
+fit(::Type{LogNormal}, ::T, ::AbstractΣstar) where T<:Real
 ```
 
 ```@docs

docs/src/partype.md

@@ -0,0 +1,66 @@
+```@meta
+CurrentModule = DistributionFits
+```
+```@meta
+DocTestSetup = :(using Statistics,DistributionFits,Optim)
+```
+
+# Creating Distributions with specific parameter types
+
+Default type of distribution parameters is Float64. Many Distributions of Distributions.jl support also other Number types, such as duals numbers or Float32 for type parameters.
+
+The fit function allow to specify the concrete parametric types of the resulting fitted distribution. Note the `Float32` parametric type in the following example.
+
+```jldoctest; output = false
+d = fit(LogNormal{Float32}, @qp_ll(1.0), @qp_uu(8))
+partype(d) == Float32
+# output
+true
+```
+
+## Infering the parametric type from other arguments.
+
+If the parametric type is omitted, default Float64 is assumed, or inferred from
+other parameters of the fitting function.
+Since quantiles and median are rather sample-like measures than paraemter-like
+measures, they do not influence the inferred parameter type.
+
+
+```jldoctest; output = false
+d = fit(LogitNormal, 0.3f0, @qp_u(0.8f0), Val(:median)) #f0 indicates Float32 literals
+partype(d) == Float64
+# output
+true
+```
+
+Parameter type, however, can be inferred from provided moments, mean, and mode.
+
+```jldoctest; output = false
+d = fit(Exponential, Moments(3f0))  #f0 indicates Float32 literals
+partype(d) == Float32
+# output
+true
+```
+
+```jldoctest; output = false
+d = fit(LogNormal, 3f0, @qp_uu(8))
+partype(d) == Float32
+# output
+true
+```
+
+```jldoctest; output = false
+d = fit(LogNormal, 3f0, @qp_uu(8), Val(:mode))
+partype(d) == Float32
+# output
+true
+```
+
+
+
+
+
+
+
+
+

src/fitstats.jl

@@ -52,6 +52,7 @@ Moments() = Moments(SA[])
 Base.getindex(m::Moments, i) = n_moments(m) >= i ? m.all[i] : 
     error("$(i)th moment not recorded.")
 Base.convert(::Type{AbstractArray}, m::Moments) = m.all
+Base.eltype(::Moments{N,T}) where {N,T} = T
 
 """
     moments(D, ::Val{N} = Val(2))
@@ -132,10 +133,10 @@ For Float64-based percentiles there are shortcuts
 - `@qp_u(q0_95)`    quantile at high p: `QuantilePoint(q0_95,0.95)`  
 - `@qp_uu(q0_975)`  quantile at very high p: `QuantilePoint(q0_975,0.975)` 
 
-For constructing QuantilePoints with type of percentiles other than Float64, 
+For constructing QuantilePoints with type of percentiles other than `Float64`, 
 use the corresponding functions, that create a percentiles of the type
-of qiven quantile.
-E.g. for a Float32-based QuantilePoint at ver low percentile
+of given quantile.
+E.g. for a `Float32`-based QuantilePoint at very low percentile
 - `qp_ll(0.2f0)` constructs a `QuantilePoint(0.2f0,0.025f0)` 
 
 There are macros/functions for some commonly used sets of QuantilePoints: 90% and 95% confidence intervals:
@@ -227,13 +228,13 @@ Fit a statistical distribution to a quantile and given statistics
 
 # Examples
 ```jldoctest fm2; output = false, setup = :(using Statistics,Distributions)
-d = fit(LogNormal, 5, @qp_uu(14));
+d = fit(LogNormal, 5.0, @qp_uu(14));
 (mean(d),quantile(d, 0.975)) .≈ (5,14)
 # output
 (true, true)
 ```
 ```jldoctest fm2; output = false, setup = :(using Statistics,Distributions)
-d = fit(LogNormal, 5, @qp_uu(14), Val(:mode));
+d = fit(LogNormal, 5.0, @qp_uu(14), Val(:mode));
 (mode(d),quantile(d, 0.975)) .≈ (5,14)
 # output
 (true, true)
@@ -245,10 +246,10 @@ function fit(::Type{D}, val, qp::QuantilePoint, ::Val{stats} = Val(:mean)) where
     stats == :median && return(fit_median_quantile(D, val, qp))
     error("unknown stats: $stats")
 end,
-function fit_median_quantile(D::Type{DT}, median, qp::QuantilePoint) where {DT <: Distribution}
+function fit_median_quantile(::Type{D}, median, qp::QuantilePoint) where {D <: Distribution}
     return(fit(D, @qp_m(median), qp))
 end,
-function fit_mean_quantile(d::Type{D}, mean::Real, qp::QuantilePoint) where D<:Distribution  
+function fit_mean_quantile(::Type{D}, mean::Real, qp::QuantilePoint) where D<:Distribution  
     error("fit_mean_quantile not yet implemented for Distribution of type: $D")
 end,
 function fit_mode_quantile(::Type{D}, mode::Real, qp::QuantilePoint)  where D<:Distribution

src/univariate/continuous/estimateMoments.jl

@@ -1,5 +1,5 @@
 # estimate the mean by numerical integration over uniform percentiles
-estimateMean(d::ContinuousUnivariateDistribution;kwargs...) = 
+estimateMean(d::ContinuousUnivariateDistribution; kwargs...) = 
   meanFunOfProb(d;kwargs...,fun=(d,p)->quantile(d,p))
 
 # estimate variance by numerical integration over uniform percentiles

src/univariate/continuous/exponential.jl

@@ -1,51 +1,69 @@
 
 
-function fit(::Type{Exponential}, m::AbstractMoments)
+fit(::Type{Exponential}, m::AbstractMoments) = fit(Exponential{eltype(m)}, m)
+function fit(::Type{Exponential{T}}, m::AbstractMoments) where T
     # https://en.wikipedia.org/wiki/Exponential_distribution
     n_moments(m) >= 1 || error("Need mean to estimate exponential")
-    return Exponential(mean(m))
+    return Exponential(T(mean(m)))
 end
 
 function fit(::Type{Exponential}, lower::QuantilePoint, upper::Missing)
+    fit(Exponential{Float64}, lower, upper)
+end
+function fit(::Type{Exponential{T}}, lower::QuantilePoint, upper::Missing) where T
     θ_lower = -lower.q/log(1-lower.p)
-    Exponential(θ_lower)
+    Exponential(T(θ_lower))
 end
 
 function fit(::Type{Exponential}, lower::Missing, upper::QuantilePoint)
+    fit(Exponential{Float64}, lower, upper)
+end
+function fit(::Type{Exponential{T}}, lower::Missing, upper::QuantilePoint) where T
     θ_upper = -upper.q/log(1-upper.p)
-    Exponential(θ_upper)
+    Exponential(T(θ_upper))
 end
 
-function fit(dt::Type{Exponential}, lower::QuantilePoint, upper::QuantilePoint)
+function fit(::Type{Exponential}, lower::QuantilePoint, upper::QuantilePoint)
+    fit(Exponential{Float64}, lower, upper)
+end
+function fit(dt::Type{Exponential{T}}, lower::QuantilePoint, upper::QuantilePoint) where T
     # return average for the two quantiles
     ismissing(lower.q) && return(fit(dt, missing, upper))
     ismissing(upper.q) && return(fit(dt, lower, missing))
     θ_lower = -lower.q/log(1-lower.p)
     θ_upper = -upper.q/log(1-upper.p)
     θ_lower ≈ θ_upper || @warn("Averaging scale for lower and upper quantile " * 
-        "for fitting expoenential distribution.")
+        "for fitting expoenential distribution.") 
     θ = (θ_lower + θ_upper)/2
-    Exponential(θ)
+    Exponential(T(θ))
 end
 
-function fit_mean_quantile(dt::Type{Exponential}, mean::Real, qp::QuantilePoint)
+function fit_mean_quantile(::Type{Exponential}, mean::T, qp::QuantilePoint) where T <: Real
+    fit_mean_quantile(Exponential{T}, mean, qp)
+end
+function fit_mean_quantile(dt::Type{Exponential{T}}, mean::Real, qp::QuantilePoint) where T
     # only fit to mean
-    warning("ignoring upper quantile when fitting Exponential to mean.")
+    @warn "ignoring upper quantile when fitting Exponential to mean." #exception=(e,catch_backtrace(),)
     fit_mean_quantile(dt, mean, missing)
 end
-
-function fit_mean_quantile(::Type{Exponential}, mean::Real, qp::Missing)
-    fit(Type{Exponential}, AbstractMoments(mean))
+function fit_mean_quantile(::Type{Exponential}, mean::T, qp::Missing) where T <: Real
+    fit(Exponential{T}, Moments(mean))
+end
+function fit_mean_quantile(::Type{Exponential{T}}, mean::Real, qp::Missing) where T <: Real
+    fit(Exponential{T}, Moments(mean))
 end
 
-function fit_mode_quantile(dt::Type{Exponential}, mode::Real, qp::QuantilePoint)
+function fit_mode_quantile(::Type{<:Exponential}, mode::T, qp::QuantilePoint) where T <: Real
     # ignore mode (its always at 0)
     mode != zero(mode) && @warn("ignoring mode when fitting Exponential.")
-    fit_mode_quantile(dt, missing, qp)
+    fit_mode_quantile(Exponential{T}, missing, qp)
 end
 
 
 function fit_mode_quantile(::Type{Exponential}, mode::Missing, qp::QuantilePoint)
+    fit_mode_quantile(Exponential{Float64}, mode, qp)
+end
+function fit_mode_quantile(::Type{Exponential{T}}, mode::Missing, qp::QuantilePoint) where T <: Real
     θ = -qp.q/log(1-qp.p)
-    Exponential(θ)
+    Exponential(T(θ))
 end

src/univariate/continuous/gamma.jl

@@ -7,6 +7,9 @@
 function fit(::Type{Gamma}, lower::QuantilePoint, upper::QuantilePoint)
     # https://www.johndcook.com/quantiles_parameters.pdf
     #    Shape-scale paras, shape α - cook:α - w:k, scale θ - cook:β - w:θ
+    if upper.p < lower.p
+        lower,upper = upper,lower
+    end
     0 < lower.p < upper.p < 1 || error(
         "Expected 0 < lower.p < upper.p < 1, but got " * 
         "lower.p = $(lower.p) and upper.p = $(upper.p)")

src/univariate/continuous/laplace.jl

@@ -1,30 +1,40 @@
 function fit(::Type{Laplace}, m::AbstractMoments)
+    fit(Laplace{eltype(m)}, m)
+end
+function fit(::Type{Laplace{T}}, m::AbstractMoments) where T
     # https://en.wikipedia.org/wiki/Laplace_distribution
     #    θ corresponds to b
     n_moments(m) >= 2 || error("Need mean and variance to estimate Laplace distribution.")
     μ = mean(m)
     θ = sqrt(var(m)/2)
-    Laplace(μ,θ)
+    Laplace(T(μ),T(θ))
 end
 
 function fit(::Type{Laplace}, lower::QuantilePoint, upper::QuantilePoint)
+    fit(Laplace{Float64}, lower, upper)
+end
+function fit(::Type{Laplace{T}}, lower::QuantilePoint, upper::QuantilePoint) where T
     # q = F^-1(p) = μ - θ a(p)
     a = (p) -> sign(p-0.5) * log(1-2*abs(p-0.5))
     a_lower, a_upper = a(lower.p), a(upper.p)
     θ = (upper.q - lower.q)/(a_lower - a_upper)
     μ = upper.q + θ * a_upper
-    Laplace(μ,θ)
+    Laplace(T(μ),T(θ))
 end
 
-function fit_mode_quantile(::Type{Laplace}, mode::Real, qp::QuantilePoint)
+function fit_mode_quantile(::Type{Laplace}, mode::T, qp::QuantilePoint) where T<:Real
+    fit_mode_quantile(Laplace{T}, mode, qp)
+end
+function fit_mode_quantile(::Type{Laplace{T}}, mode::Real, qp::QuantilePoint) where T
     a = (p) -> sign(p-0.5) * log(1-2*abs(p-0.5))
     a_qp = a(qp.p)
     #θ = (qp.q - mode)/(0 - a_qp)
     θ = (mode - qp.q)/a_qp
-    Laplace(mode,θ)
+    Laplace(T(mode),T(θ))
 end
 
-function fit_mean_quantile(::Type{Laplace}, mean::Real, qp::QuantilePoint)
-    fit_mode_quantile(Laplace, mean, qp)
+function fit_mean_quantile(D::Type{<:Laplace}, mean::Real, qp::QuantilePoint)
+    # mode equals mean
+    fit_mode_quantile(D, mean, qp)
 end