diff --git a/docs/src/api.md b/docs/src/api.md index 7b6a272e..aa8f17e4 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -46,6 +46,8 @@ atan sinh cosh tanh +gaussian +erf abs ``` @@ -72,5 +74,7 @@ acosHomogCoef atanHomogCoef sinhcoshHomogCoef tanhHomogCoef +gaussHomogCoef +erfHomogCoef A_mul_B! ``` diff --git a/src/Taylor1.jl b/src/Taylor1.jl index fcafb423..9d8adef5 100644 --- a/src/Taylor1.jl +++ b/src/Taylor1.jl @@ -1118,6 +1118,138 @@ function tanhHomogCoef{T<:Number}(kcoef::Int,ac::Array{T,1},coeffst2::Array{T,1} coefhomog end + +### STAT & MISC FUNCTIONS ### + +## gaussian ## +doc""" + gaussian(a) + +Return the Taylor expansion of $e^{-a^2}$, of order `a.order`, for +`a::Taylor1` polynomial. + +For details on making the Taylor expansion, see [`TaylorSeries.gaussHomogCoef`](@ref). +""" +function gaussian(a::Taylor1) + @inbounds aux = exp(-a.coeffs[1]^2) + T = typeof(aux) + a_prime = derivative(a) + ac_prime = convert(Vector{T},a_prime.coeffs) + ac = convert(Vector{T},a.coeffs) + coeffs = zeros(ac) + kcoefprod = zeros(ac) + @inbounds coeffs[1] = convert(T,aux) + @inbounds kcoefprod[1] = convert(T,mulHomogCoef(0,ac_prime,ac)) + @inbounds for k in 1:a.order-1 + kcoefprod[k+1] = mulHomogCoef(k,ac_prime,ac) + coeffs[k+1] = gaussHomogCoef(k,kcoefprod,coeffs) + end + Taylor1( coeffs, a.order) +end + +## standarized gaussian ## +doc""" + gaussian(μ,σ,order) + +Return the Taylor expansion of $\frac{1}{\sigma \sqrt{2 \pi} } e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^2} $, of order `order`, for +`a::Taylor1` polynomial, where μ is the median and σ the standard deviation. + +For details on making the Taylor expansion, see [`TaylorSeries.gaussHomogCoef`](@ref). +""" +function gaussian(μ=0.,σ=1./sqrt(2),order=15) + a = Taylor1(order) + a = (a-μ)/(sqrt(2.)*σ) + @inbounds aux = exp(-a.coeffs[1]^2) + T = typeof(aux) + a_prime = derivative(a) + ac_prime = convert(Vector{T},a_prime.coeffs) + ac = convert(Vector{T},a.coeffs) + coeffs = zeros(ac) + kcoefprod = zeros(ac) + @inbounds coeffs[1] = convert(T,aux) + @inbounds kcoefprod[1] = convert(T,mulHomogCoef(0,ac_prime,ac)) + @inbounds for k in 1:a.order-1 + kcoefprod[k+1] = mulHomogCoef(k,ac_prime,ac) + coeffs[k+1] = gaussHomogCoef(k,kcoefprod,coeffs) + end + 1./(σ*sqrt(2*pi)) * Taylor1( coeffs, a.order) +end + +# Homogeneous coefficients for gaussian +doc""" +gaussHomogCoef(kcoef, kcoefprod, coeffs) + +Compute the `k-th` expansion coefficient of $g(a) = e^{-a^2}$ given by + +\begin{equation*} +g_k = -\frac{2}{k} \sum_{j=0}^{k-1} p_{k-j} g_j, +\end{equation*} + +with $a$ a `Taylor1` polynomial and $ p(x) = a(x) a'(x) $. + +Inputs are the `kcoef`-th coefficient, the already calculated expansion coefficients `kcoefprod` of $p(x)$ and the already calculated expansion coefficients `coeffs` of $g(a)$. +""" +function gaussHomogCoef{T<:Number}(kcoef::Int,kcoefprod::Vector{T},coeffs::Vector{T}) + #kcoef == 0 && return exp(-coeffs[1]^2) + coefhomog = zero(T) + @inbounds for j = 0:kcoef-1 + coefhomog += kcoefprod[j+1] * coeffs[kcoef-j] + end + -2*coefhomog / (kcoef) +end + +## erf ## +doc""" + erf(a) + +Return the Taylor expansion of $\frac{2}{\sqrt{\pi}}\int_0^a{e^{-t^2}dt}$, of order `a.order`, for +`a::Taylor1` polynomial. + +For details on making the Taylor expansion, see [`TaylorSeries.gaussHomogCoef`](@ref). +""" +function erf(a::Taylor1) + @inbounds aux = sqrt(pi)/2 * erf(a.coeffs[1]) + T = typeof(aux) + a_prime = derivative(a) + ac_prime = convert(Vector{T},a_prime.coeffs) + ac = convert(Vector{T},a.coeffs) + coeffs = zeros(ac) + kcoefprod = zeros(ac) + gausscoeffs = zeros(ac) + @inbounds coeffs[1] = convert(T,aux) + @inbounds kcoefprod[1] = convert(T,mulHomogCoef(0,ac_prime,ac)) + @inbounds gausscoeffs[1] = convert(T,exp(-a.coeffs[1]^2)) + @inbounds for k in 1:a.order + kcoefprod[k+1] = mulHomogCoef(k,ac_prime,ac) + gausscoeffs[k+1] = gaussHomogCoef(k,kcoefprod,gausscoeffs) + coeffs[k+1] = erfHomogCoef(k,gausscoeffs,ac_prime) + end + 2/sqrt(pi) * Taylor1( coeffs, a.order) +end + +# Homogeneous coefficients for erf +doc""" + erfHomogCoef(kcoef, gausscoeffs, coeffs) + +Compute the `k-th` expansion coefficient of $e(a) = erf(a)$ given by + +\begin{equation*} +e_k = \frac{1}{k} \sum_{j=0}^{k-1} g_{k-j} a'_j, +\end{equation*} + +with $a$ a `Taylor1` polynomial, $g(a) = e^{-a^2}$ and $a'(x) = \frac{\partial{}}{\partial{x}}a$. + +Inputs are the `kcoef`-th coefficient, the already calculated expansion coefficients `gausscoeffs` of $g(a)$ and the already calculated expansion coefficients `coeffs` of $e(a)$. +""" +function erfHomogCoef{T<:Number}(kcoef::Int,gausscoeffs::Vector{T},primecoeffs::Vector{T}) + #kcoef == 0 && return exp(-coeffs[1]^2) + coefhomog = zero(T) + @inbounds for j = 0:kcoef-1 + coefhomog += primecoeffs[j+1] * gausscoeffs[kcoef-j] + end + coefhomog / (kcoef) +end + ## Differentiating ## """ derivative(a) diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl index 00db7005..5c2996ae 100644 --- a/src/TaylorSeries.jl +++ b/src/TaylorSeries.jl @@ -26,11 +26,12 @@ import Base: zero, one, zeros, ones, isinf, isnan, rem, mod, mod2pi, abs, gradient, sqrt, exp, log, sin, cos, tan, asin, acos, atan, sinh, cosh, tanh, + erf, reverse, A_mul_B! export Taylor1, TaylorN, HomogeneousPolynomial -export get_coeff, derivative, integrate, +export get_coeff, derivative, gaussian, integrate, evaluate, evaluate!, show_params_TaylorN, get_order, get_numvars, diff --git a/test/runtests.jl b/test/runtests.jl index 454669d9..b315176b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -127,6 +127,12 @@ facts("Tests for Taylor1 expansions") do @fact evaluate(tanh(t/2),1.5) == evaluate(sinh(t) / (cosh(t) + 1),1.5) --> true @fact cosh(t) == real(cos(im*t)) --> true @fact sinh(t) == imag(sin(im*t)) --> true + + μ , σ = (3.2, 1) + @fact gaussian(μ,σ,t.order) == 1/(sqrt(2pi)*σ)*gaussian((t-μ)/(sqrt(2.)*σ)) --> true + @fact gaussian(t) == exp(-t^2) --> true + @fact erf(t) == 2/sqrt(pi) * integrate(exp(-t^2)) --> true + @fact erf(t) == 2/sqrt(pi) * integrate( gaussian(t) ) --> true v = [sin(t), exp(-t)] @fact evaluate(v) == [0.0, 1.0] --> true