diff --git a/docs/src/fit.md b/docs/src/fit.md
index 24e8ee1dc..eb91aa39b 100644
--- a/docs/src/fit.md
+++ b/docs/src/fit.md
@@ -53,6 +53,7 @@ The `fit_mle` method has been implemented for the following distributions:
 - [`Rayleigh`](@ref)
 - [`InverseGaussian`](@ref)
 - [`Uniform`](@ref)
+- [`Weibull`](@ref)
 
 **Multivariate:**
 
diff --git a/src/univariate/continuous/weibull.jl b/src/univariate/continuous/weibull.jl
index 40f50533a..d081c0717 100644
--- a/src/univariate/continuous/weibull.jl
+++ b/src/univariate/continuous/weibull.jl
@@ -139,3 +139,55 @@ end
 #### Sampling
 
 rand(rng::AbstractRNG, d::Weibull) = xv(d, randexp(rng))
+
+#### Fit model
+
+"""
+    fit_mle(::Type{<:Weibull}, x::AbstractArray{<:Real}; 
+    alpha0::Real = 1, maxiter::Int = 1000, tol::Real = 1e-16)
+
+Compute the maximum likelihood estimate of the [`Weibull`](@ref) distribution with Newton's method.
+"""
+function fit_mle(::Type{<:Weibull}, x::AbstractArray{<:Real};
+    alpha0::Real = 1, maxiter::Int = 1000, tol::Real = 1e-16)
+
+    N = 0
+
+    lnx = map(log, x)
+    lnxsq = lnx.^2
+    mean_lnx = mean(lnx)
+
+    # first iteration outside loop, prevents type instabililty in α, ϵ
+
+    xpow0 = x.^alpha0
+    sum_xpow0 = sum(xpow0)
+    dot_xpowlnx0 = dot(xpow0, lnx)
+
+    fx = dot_xpowlnx0 / sum_xpow0 - mean_lnx - 1 / alpha0
+    ∂fx = (-dot_xpowlnx0^2 + sum_xpow0 * dot(lnxsq, xpow0)) / (sum_xpow0^2) + 1 / alpha0^2
+
+    Δα = fx / ∂fx
+    α = alpha0 - Δα
+
+    ϵ = abs(Δα)
+    N += 1
+
+    while ϵ > tol && N < maxiter
+
+        xpow = x.^α
+        sum_xpow = sum(xpow)
+        dot_xpowlnx = dot(xpow, lnx)
+
+        fx = dot_xpowlnx / sum_xpow - mean_lnx - 1 / α
+        ∂fx = (-dot_xpowlnx^2 + sum_xpow * dot(lnxsq, xpow)) / (sum_xpow^2) + 1 / α^2
+
+        Δα = fx / ∂fx
+        α -= Δα
+
+        ϵ = abs(Δα)
+        N += 1
+    end
+
+    θ = mean(x.^α)^(1 / α)
+    return Weibull(α, θ)
+end
diff --git a/test/fit.jl b/test/fit.jl
index 4cf299c6a..f774b4191 100644
--- a/test/fit.jl
+++ b/test/fit.jl
@@ -443,3 +443,13 @@ end
         @test all(ForwardDiff.gradient(f, x) .>= 0)
     end
 end
+
+@testset "Testing fit for Weibull" begin
+    for func in funcs, dist in (Weibull, Weibull{Float64})
+        d = fit(dist, func[2](dist(8.1, 4.3), N))
+        @test isa(d, dist)
+        @test isapprox(d.α, 8.1, atol = 0.1)
+        @test isapprox(d.θ, 4.3, atol = 0.1)
+
+    end
+end