As a suggestion, if you want to try to leverage all of Julia, I think something like Reduce is the way to go.

For example, if you want to calculate partial Normal moments

```
julia> using Reduce, Compat, SpecialFunctions
Reduce (Free PSL version, revision 4218), 22-Sep-2017 ...
julia> Reduce.Rational(false);
julia> @generated function normal_moment(x, ::Val{N} = Val{0}()) where N
ex = int(:( x^$N * exp(-x^2/2 ) ), :x)
root_2π = sqrt(2π)
:( $ex / $root_2π )
end
julia> normal_moment(0)
0.0
julia> normal_moment(0.5)
0.19146246127401312
julia> normal_moment(0.5, Val(2))
0.015429797891863349
julia> normal_moment(1.5, Val(2))
0.23891640523230437
julia> normal_moment(10, Val(4)) ###Pick a big number to get the kurtosis of a standard normal.
1.5000000000000002
```

This could be useful if you wanted to integrate Hermite polynomials.

If you like statistics, you can also have a lot of fun with moment generating functions:

```
julia> @generated function gamma_moment(a, b, ::Val{N} = Val{0}()) where N
ex = :((1 - t/b)^-a)
for i ∈ 1:N
ex = df( ex, :t )
end
quote
t = 0
$ex
end
end
gamma_moment (generic function with 2 methods)
julia> gamma_moment(2, 2, Val(1)) ##The mean
1.0
julia> gamma_var(α, β) = gamma_moment(α, β, Val(2)) - gamma_moment(α, β, Val(1))
1.5
julia> gamma_var(2, 2)
0.5
julia> gamma_var(3.4, 1.7)
1.1764705882352953
julia> gamma_var(4, 2)
1.0
```

Be warned though that if the Val you’re using isn’t known at compile time, the code wont be type stable.