diff --git a/Project.toml b/Project.toml
index 643ea0f47..e4bf95ef4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,6 +9,8 @@ ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 DensityInterface = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6"
 PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
@@ -43,6 +45,8 @@ FiniteDifferences = "0.12"
 ForwardDiff = "0.10, 1"
 JSON = "0.21"
 LinearAlgebra = "<0.0.1, 1"
+LogExpFunctions = "0.3.29"
+IrrationalConstants = "0.2.6"
 OffsetArrays = "1"
 PDMats = "0.11.35"
 Printf = "<0.0.1, 1"
diff --git a/src/Distributions.jl b/src/Distributions.jl
index fcf699434..cfe5204fc 100644
--- a/src/Distributions.jl
+++ b/src/Distributions.jl
@@ -25,8 +25,9 @@ import StatsBase: kurtosis, skewness, entropy, mode, modes,
 import PDMats: dim, PDMat, invquad
 
 using SpecialFunctions
+using LogExpFunctions: logabssinh, log1psq
 using Base.MathConstants: eulergamma
-
+using IrrationalConstants: invπ
 import AliasTables
 
 export
diff --git a/src/univariate/continuous/cosine.jl b/src/univariate/continuous/cosine.jl
index 4d11a9ddd..ae068ccac 100644
--- a/src/univariate/continuous/cosine.jl
+++ b/src/univariate/continuous/cosine.jl
@@ -68,7 +68,14 @@ function pdf(d::Cosine{T}, x::Real) where T<:Real
     end
 end
 
-logpdf(d::Cosine, x::Real) = log(pdf(d, x))
+function logpdf(d::Cosine{T}, x::Real) where T<:Real
+    if insupport(d, x)
+        z = (x - d.μ) / d.σ
+        return log1p(cospi(z)) - log(2d.σ)
+    else
+        return typemin(T)
+    end
+end
 
 function cdf(d::Cosine{T}, x::Real) where T<:Real
     if x < d.μ - d.σ
@@ -93,3 +100,27 @@ function ccdf(d::Cosine{T}, x::Real) where T<:Real
 end
 
 quantile(d::Cosine, p::Real) = quantile_bisect(d, p)
+
+function mgf(d::Cosine, t::Real)
+    σt, μt = d.σ * t, d.μ * t
+    z = iszero(σt) ? one(float(σt)) : sinh(σt)/σt
+    return exp(μt) * (z / (1 + (invπ * σt)^2))
+end
+
+function cgf(d::Cosine, t::Real)
+    σt, μt = d.σ * t, d.μ * t
+    z = iszero(σt) ? zero(float(σt)) : logabssinh(σt) - log(σt)
+    return μt + z - log1psq(invπ * σt)
+end
+
+function cf(d::Cosine, t::Real)
+    σt, μt = d.σ * t, d.μ * t
+    abs(σt) ≈ π && return cis(μt) / 2
+    z = iszero(σt) ? one(float(σt)) : sin(σt)/σt
+    return cis(μt) * z / (1 - (invπ * σt)^2)
+end
+
+#### Affine transformations
+
+Base.:+(d::Cosine, a::Real) = Cosine(d.μ + a, d.σ)
+Base.:*(c::Real, d::Cosine) = Cosine(c * d.μ, abs(c) * d.σ)
diff --git a/test/univariate/continuous/cosine.jl b/test/univariate/continuous/cosine.jl
new file mode 100644
index 000000000..5872cf3a1
--- /dev/null
+++ b/test/univariate/continuous/cosine.jl
@@ -0,0 +1,14 @@
+@testset "Cosine" begin
+    @testset "mgf, cgf, cf" begin
+        @test cf(Cosine(1, 1), 1) ≈ cis(1) * sin(1)/(1 - π^-2)
+        @test mgf(Cosine(-1, 2), -1) ≈ exp(1) * sinh(-2) / (-2*(1 + 4π^-2))
+        test_cgf(Cosine(-1, 2), [-1,2, 4])
+    end
+
+    @testset "affine transformations" begin
+        test_affine_transformations(Cosine, (4, 2))
+        @test -Cosine(1, 1) == Cosine(-1, 1)
+        @test 3Cosine(1, 2) == Cosine(3, 6)
+        @test Cosine(0, 2) + 42 == Cosine(42, 2)
+    end
+end