diff --git a/Project.toml b/Project.toml
index 85684e085c..5569a6b724 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,13 +1,14 @@
 name = "Distributions"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
-version = "0.25.123"
 authors = ["JuliaStats"]
+version = "0.25.124"
 
 [deps]
 AliasTables = "66dad0bd-aa9a-41b7-9441-69ab47430ed8"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 DensityInterface = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
@@ -42,6 +43,7 @@ Distributed = "<0.0.1, 1"
 FillArrays = "0.9, 0.10, 0.11, 0.12, 0.13, 1"
 FiniteDifferences = "0.12"
 ForwardDiff = "0.10, 1"
+HypergeometricFunctions = "0.3.28"
 JSON = "0.21, 1"
 LinearAlgebra = "<0.0.1, 1"
 OffsetArrays = "1"
diff --git a/src/Distributions.jl b/src/Distributions.jl
index b18dfa8086..c7e453762d 100644
--- a/src/Distributions.jl
+++ b/src/Distributions.jl
@@ -25,6 +25,7 @@ import StatsBase: kurtosis, skewness, entropy, mode, modes,
 import PDMats: dim, PDMat, invquad
 
 using SpecialFunctions
+import HypergeometricFunctions
 using Base.MathConstants: eulergamma
 
 import AliasTables
diff --git a/src/univariate/continuous/rician.jl b/src/univariate/continuous/rician.jl
index f5c354f705..0cf842be06 100644
--- a/src/univariate/continuous/rician.jl
+++ b/src/univariate/continuous/rician.jl
@@ -63,11 +63,11 @@ partype(d::Rician{T}) where {T<:Real} = T
 
 #### Statistics
 
-# helper
-_Lhalf(x) = exp(x/2) * ((1-x) * besseli(zero(x), -x/2) - x * besseli(oneunit(x), -x/2))
+# Laguerre function L_{1/2}
+_Lhalf(x::Real) = HypergeometricFunctions._₁F₁(-one(x)/2, one(x), x)
 
-mean(d::Rician) = d.σ * sqrthalfπ * _Lhalf(-d.ν^2/(2 * d.σ^2))
-var(d::Rician) = 2 * d.σ^2 + d.ν^2 - halfπ * d.σ^2 * _Lhalf(-d.ν^2/(2 * d.σ^2))^2
+mean(d::Rician) = d.σ * sqrthalfπ * _Lhalf(-(d.ν / d.σ)^2/2)
+var(d::Rician) = 2 * d.σ^2 + d.ν^2 - halfπ * d.σ^2 * _Lhalf(-(d.ν/d.σ)^2/2)^2
 
 function mode(d::Rician)
     m = mean(d)
@@ -91,7 +91,21 @@ function _minimize_gss(f, a, b; tol=1e-12)
     (b + a) / 2
 end
 
-#### PDF/CDF/quantile delegated to NoncentralChisq
+#### PDF
+
+function pdf(d::Rician, x::Real)
+    ν, σ = params(d)
+    result = x / σ^2 * exp(-((x-ν)/σ)^2/2) * besselix(0, (x * ν)/ σ^2)
+    return x < 0 || isinf(x) ? zero(result) : result
+end
+
+function logpdf(d::Rician, x::Real)
+    ν, σ = params(d)
+    result = log(abs(x) / σ^2 * besselix(0, (x * ν)/ σ^2)) - ((x-ν)/σ)^2/2
+    return x < 0 || isinf(x) ? oftype(result, -Inf) : result
+end
+
+#### quantile/CDF delegated to NoncentralChisq
 
 function quantile(d::Rician, x::Real)
     ν, σ = params(d)
@@ -103,16 +117,14 @@ function cquantile(d::Rician, x::Real)
     return sqrt(cquantile(NoncentralChisq(2, (ν / σ)^2), x)) * σ
 end
 
-function pdf(d::Rician, x::Real)
+function invlogcdf(d::Rician, lx::Real)
     ν, σ = params(d)
-    result = 2 * x / σ^2 * pdf(NoncentralChisq(2, (ν / σ)^2), (x / σ)^2)
-    return x < 0 || isinf(x) ? zero(result) : result
+    return sqrt(invlogcdf(NoncentralChisq(2, (ν / σ)^2), lx)) * σ
 end
 
-function logpdf(d::Rician, x::Real)
+function invlogccdf(d::Rician, lx::Real)
     ν, σ = params(d)
-    result = log(2 * abs(x) / σ^2) + logpdf(NoncentralChisq(2, (ν / σ)^2), (x / σ)^2)
-    return x < 0 || isinf(x) ? oftype(result, -Inf) : result
+    return sqrt(invlogccdf(NoncentralChisq(2, (ν / σ)^2), lx)) * σ
 end
 
 function cdf(d::Rician, x::Real)
@@ -127,6 +139,18 @@ function logcdf(d::Rician, x::Real)
     return x < 0 ? oftype(result, -Inf) : result
 end
 
+function ccdf(d::Rician, x::Real)
+    ν, σ = params(d)
+    result = ccdf(NoncentralChisq(2, (ν / σ)^2), (x / σ)^2)
+    return x < 0 ? one(result) : result
+end
+
+function logccdf(d::Rician, x::Real)
+    ν, σ = params(d)
+    result = logccdf(NoncentralChisq(2, (ν / σ)^2), (x / σ)^2)
+    return x < 0 ? zero(result) : result
+end
+
 #### Sampling
 
 function rand(rng::AbstractRNG, d::Rician)
diff --git a/test/univariate/continuous/rician.jl b/test/univariate/continuous/rician.jl
index a75397f891..5b5bc886a3 100644
--- a/test/univariate/continuous/rician.jl
+++ b/test/univariate/continuous/rician.jl
@@ -1,3 +1,5 @@
+using StatsFuns
+
 @testset "Rician" begin
 
     d1 = Rician(0.0, 10.0)
@@ -9,6 +11,24 @@
     @test eltype(d1) === Float64
     @test rand(d1) isa Float64
 
+    # Reference: WolframAlpha
+    @test @inferred(mean(d1))::Float64 ≈ 5 * sqrt(2 * π)
+    @test @inferred(std(d1))::Float64 ≈ sqrt(200 - 50 * π)
+    @test @inferred(var(d1))::Float64 ≈ 200 - 50 * π
+    @test @inferred(mode(d1))::Float64 ≈ 10
+    ps = 0:0.1:1
+    @test all(splat(isapprox), zip(quantile.(d1, ps), @. sqrt(-200 * log1p(-ps))))
+    @test all(splat(isapprox), zip(cquantile.(d1, ps), @. sqrt(-200 * log(ps))))
+    @test all(splat(isapprox), zip(invlogcdf.(d1, log.(ps)), @. sqrt(-200 * log1p(-ps))))
+    @test all(splat(isapprox), zip(invlogccdf.(d1, log.(ps)), @. sqrt(-200 * log(ps))))
+    xs = 0:0.5:30 # 99th percentile is approx 30
+    @test all(splat(isapprox), zip(pdf.(d1, xs), map(x -> x / 100 * exp(-x^2/200), xs)))
+    @test all(splat(isapprox), zip(logpdf.(d1, xs), map(x -> log(x / 100) - x^2/200, xs)))
+    @test all(splat(isapprox), zip(cdf.(d1, xs), @. 1 - exp(-xs^2/200)))
+    @test all(splat(isapprox), zip(logcdf.(d1, xs), @. log1mexp(-xs^2/200)))
+    @test all(splat(isapprox), zip(ccdf.(d1, xs), @. exp(-xs^2/200)))
+    @test all(splat(isapprox), zip(logccdf.(d1, xs), @. -xs^2/200))
+
     d2 = Rayleigh(10.0)
     @test mean(d1) ≈ mean(d2)
     @test var(d1) ≈ var(d2)
@@ -23,7 +43,22 @@
     x = Base.Fix1(quantile, d1).([0.25, 0.45, 0.60, 0.80, 0.90])
     @test all(Base.Fix1(cdf, d1).(x) .≈ [0.25, 0.45, 0.60, 0.80, 0.90])
 
-    x = rand(Rician(5.0, 5.0), 100000)
+    # Reference: WolframAlpha
+    @test mean(d1) ≈ 15.4857246055
+    @test std(d1) ≈ 7.75837182934
+    @test var(d1) ≈ 60.1923334423
+    @test median(d1) ≈ 14.7547909179
+    @test all(xy -> isapprox(xy...; rtol=1e-5), zip(quantile.(d1, [0.1, 0.25, 0.5, 0.75, 0.9]), [5.86808, 9.62923, 14.7548, 20.5189, 26.0195]))
+
+    d1 = Rician(5.0, 5.0)
+    # Reference: WolframAlpha
+    @test mean(d1) ≈ 7.7428623027557
+    @test std(d1) ≈ 3.8791859146687
+    @test var(d1) ≈ 15.0480833605606
+    @test median(d1) ≈ 7.37739545894
+    @test all(xy -> isapprox(xy...; rtol=1e-5), zip(quantile.(d1, [0.1, 0.25, 0.5, 0.75, 0.9]), [2.93404, 4.81462, 7.3774, 10.2595, 13.0097]))
+
+    x = rand(d1, 100000)
     d1 = fit(Rician, x)
     @test d1 isa Rician{Float64}
     @test params(d1)[1] ≈ 5.0 atol=0.2
@@ -37,6 +72,11 @@
     @test partype(d1) === Float32
     @test eltype(d1) === Float64
     @test rand(d1) isa Float64
+    d1_f64 = Rician(10.0, 10.0)
+    @test @inferred(mean(d1))::Float32 ≈ Float32(mean(d1_f64))
+    @test @inferred(std(d1))::Float32 ≈ Float32(std(d1_f64))
+    @test @inferred(var(d1))::Float32 ≈ Float32(var(d1_f64))
+    @test @inferred(median(d1))::Float32 ≈ Float32(median(d1_f64))
 
     d1 = Rician()
     @test d1 isa Rician{Float64}
@@ -94,4 +134,19 @@
     @inferred logcdf(d1, 1//2)
     @inferred logcdf(d1, Inf)
 
+    # issue #2035
+    # reference values computed with WolframAlpha
+    d = Rician(1.0, 0.12)
+    @test mean(d) ≈ 1.00722650704
+    @test std(d) ≈ 0.119560710568
+    @test var(d) ≈ 0.0142947635114
+    @test median(d) ≈ 1.00719144719
+    @test pdf(d, 0.5) ≈ 0.000400758853946
+    @test pdf(d, 1) ≈ 3.33055235263
+    @test pdf(d, 1.5) ≈ 0.000692437774661
+    @test pdf(d, 2) ≈ 3.91711741719e-15
+    @test logpdf(d, 0.5) ≈ -7.82215067328
+    @test logpdf(d, 1) ≈ 1.2031381619
+    @test logpdf(d, 1.5) ≈ -7.27529218003
+    @test logpdf(d, 2) ≈ -33.1734203644
 end