Fix fit docs (#1746)

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
Co-authored-by: Viral B. Shah <ViralBShah@users.noreply.github.com>

docs/src/fit.md

@@ -27,6 +27,8 @@ The function `fit_mle` is for maximum likelihood estimation.
 ### Synopsis
 
 ```@docs
+fit(D, x)
+fit(D, x, w)
 fit_mle(D, x)
 fit_mle(D, x, w)
 ```

src/genericfit.jl

@@ -30,5 +30,18 @@ fit_mle(dt::Type{D}, x::AbstractArray, w::AbstractArray) where {D<:UnivariateDis
 fit_mle(dt::Type{D}, x::AbstractMatrix) where {D<:MultivariateDistribution} = fit_mle(D, suffstats(D, x))
 fit_mle(dt::Type{D}, x::AbstractMatrix, w::AbstractArray) where {D<:MultivariateDistribution} = fit_mle(D, suffstats(D, x, w))
 
+"""
+    fit(D, args...)
+
+Fit a distribution of type `D` to `args`.
+
+The fit function will choose a reasonable way to fit the distribution, which,
+in most cases, is maximum likelihood estimation. Note that this algorithm may
+change; for a function that will behave consistently across versions, see 
+`fit_mle`.
+
+By default, the fallback is [`fit_mle(D, args...)`](@ref); developers can change this default
+for a specific distribution type `D <: Distribution` by defining a `fit(::Type{D}, args...)` method.
+"""
 fit(dt::Type{D}, x) where {D<:Distribution} = fit_mle(D, x)
 fit(dt::Type{D}, args...) where {D<:Distribution} = fit_mle(D, args...)

src/univariate/continuous/beta.jl

@@ -208,12 +208,8 @@ function rand(rng::AbstractRNG, d::Beta{T}) where T
     end
 end
 
-#### Fit model
-"""
-    fit_mle(::Type{<:Beta}, x::AbstractArray{T})
 
-Maximum Likelihood Estimate of `Beta` Distribution via Newton's Method
-"""
+
 function fit_mle(::Type{<:Beta}, x::AbstractArray{T};
     maxiter::Int=1000, tol::Float64=1e-14) where T<:Real
 
@@ -240,11 +236,8 @@ function fit_mle(::Type{<:Beta}, x::AbstractArray{T};
     return Beta(θ[1], θ[2])
 end
 
-"""
-    fit(::Type{<:Beta}, x::AbstractArray{T})
 
-fit a `Beta` distribution
-"""
+
 function fit(::Type{<:Beta}, x::AbstractArray{T}) where T<:Real
     x_bar = mean(x)
     v_bar = varm(x, x_bar)