function gauss(p::Vector{Float64}; x::Float64=0.0)
  λ, μ, σ = p
  return λ / sqrt(2π * σ^2) * exp(-(x - μ)^2 / 2σ^2)
end

function J_gauss(p::Vector{Float64}; x::Float64=0.0)
  λ, μ, σ = p
  c = [(1 / λ) ((x - μ) / σ^2) ((x - μ)^2 - σ^2) / σ^3]
  return c .* gauss(p; x=x)
end

function H_gauss(p::Vector{Float64}; x::Float64=0.0)
  λ, μ, σ = p
  c = zeros(3, 3)

  # c[1,1] = 0
  c[1, 2] = c[2, 1] = (x - μ) / (λ * σ^2)
  c[1, 3] = c[3, 1] = -((μ + σ - x) * (-μ + σ + x)) / (λ * σ^3)

  c[2, 2] = ((x - μ)^2 - σ^2) / σ^4
  c[2, 3] = c[3, 2] = -((μ - x) * ((μ - x)^2 - 3σ^2)) / σ^5

  c[3, 3] = (-5 * σ^2 * (μ - x)^2 + (μ - x)^4 + 2 * σ^4) / σ^6
  return c .* gauss(p; x=x)
end

@testset "Gauss function" begin
  λ = 10 * rand()
  μ = 10 * (rand() - 0.5)
  σ = 10rand()
  p = [λ, μ, σ]

  x = σ * (2 * rand() - 1) + μ

  @assert abs(gauss(p; x=x)) > 1e-7 "gauss(p; x=x) is too small: $(gauss(p; x=x))"

  @test TaylorTest.check(gauss, J_gauss, p; x=x)
  @test TaylorTest.check(J_gauss, H_gauss, p; x=x)
end

@testset "Composed Gauss function" begin
  p = 1 .+ 5rand(6)

  p_aux = _p -> [_p[1] / _p[2], _p[3] + _p[4] + 2_p[5], _p[6]]

  J_p_aux = _p -> [1/_p[2] -_p[1]/_p[2]^2 0 0 0 0;
    0 0 1 1 2 0;
    0 0 0 0 0 1]

  H_p_aux = _p -> begin
    h = zeros(3, 6, 6)
    h[1, 1, 2] = h[1, 2, 1] = -1 / _p[2]^2
    h[1, 2, 2] = 2_p[1] / _p[2]^3
    h
  end

  g = (_p; x = 0.0) -> gauss(p_aux(_p); x=x)
  Jg = (_p; x = 0.0) -> J_gauss(p_aux(_p); x=x) * J_p_aux(_p)
  Hg = (_p; x = 0.0) -> begin
    p = p_aux(_p)
    _J_gauss = vec(J_gauss(p; x=x))
    _H_p = H_p_aux(_p)
    _J_p = J_p_aux(_p)
    r = _J_p' * H_gauss(p; x=x) * _J_p
    r .+= TO.tensorcontract([2, 3], _J_gauss, [1], _H_p, [1, 2, 3])
    r
  end

  _p = p_aux(p)
  x = _p[2] + _p[2] * 2 * (rand() - 0.5)

  @test TaylorTest.check(g, Jg, p; x=x)
  @test TaylorTest.check(Jg, Hg, p; x=x)
end