Problem in using turing_inference() for Lorenz equation.

[Vaibhavdixit02]

When using turing_inference for Lorenz equation the priors seem to be causing some problems, in case of a model like

g1 = @ode_def_bare LorenzExample begin
  dx = σ*(y-x)
  dy = x*(ρ-z) - y
  dz = x*y - β*z
end σ ρ β 
r0 = [1.0; 0.0; 0.0]                
tspan = (0.0, 30.0)
p = [10.0,28.0,2.66]
prob = ODEProblem(g1,r0,tspan,p)
@time sol = solve(prob,Vern9(),abstol=1e-12,reltol=1e-12)

t = collect(linspace(1,30,30))
sig = 0.49
data = convert(Array, VectorOfArray([(sol(t[i]) + sig*randn(3)) for i in 1:length(t)]))

priors = [Truncated(Normal(10,2),0,15),Truncated(Normal(30,5),0,45),Truncated(Normal(2.5,0.5),0,4)]

The turing_inference call

@time bayesian_result = turing_inference(prob,Tsit5(),t,data,priors;num_samples=500)

Gives the following error repeatedly

[Turing.WARNING]: Numerical error has been found in gradients.
 verifygrad(::Array{Float64,1}) at ad.jl:87

@xukai92 can you shed some light on why this is happening? Not using Truncated priors does not give this error but the result is really bad in that case.

[Vaibhavdixit02]

@xukai92 could you please take a look at this. Thanks!

[ChrisRackauckas]

Seems like it could be HMC errors. Can you try manually transforming the parameters to use an exponential? Is there a way to set this in Turing.jl?

[xukai92]

I think it should be the step size of HMC is set to a too large value.

[xukai92]

Emmm but your priors seem to be fine.

[xukai92]

Using a smaller step size, i.e. @time bayesian_result = turing_inference(prob,Tsit5(),t,data,priors;num_samples=500,epsilon = 0.001) works for me

[ChrisRackauckas]

When using that it still throws

[Turing.WARNING]: Numerical error has been found in gradients.
 verifygrad(::Array{Float64,1}) at ad.jl:87

for me, though it runs (but gets the incorrect result).

[xukai92]

It might caused by a bug I solved in TuringLang/Turing.jl@170c2af about initialization.

BTW how should I check the results? I can play with it to make sure it works.

[ChrisRackauckas]

using DiffEqBayes, OrdinaryDiffEq, ParameterizedFunctions, RecursiveArrayTools
g1 = @ode_def_bare LorenzExample begin
  dx = σ*(y-x)
  dy = x*(ρ-z) - y
  dz = x*y - β*z
end σ ρ β
r0 = [1.0; 0.0; 0.0]
tspan = (0.0, 30.0)
p = [10.0,28.0,2.66]
prob = ODEProblem(g1,r0,tspan,p)
@time sol = solve(prob,Vern9(),abstol=1e-12,reltol=1e-12)

t = collect(linspace(1,30,30))
sig = 0.49
data = convert(Array, VectorOfArray([(sol(t[i]) + sig*randn(3)) for i in 1:length(t)]))

priors = [Truncated(Normal(10,2),0,15),Truncated(Normal(30,5),0,45),Truncated(Normal(2.5,0.5),0,4)]

@time bayesian_result = turing_inference(prob,Tsit5(),t,data,priors;num_samples=500,epsilon = 0.001)

works well for me now. I don't know why it's different from before. @Vaibhavdixit02 try it?

[Vaibhavdixit02]

The warning still appears for me despite multiple efforts, @xukai92 can you try with a version of Turing.jl before the mentioned PR?

[xukai92]

@Vaibhavdixit02 I tried it with current release Turing.jl and it works with randomness. Just do multiple runs will always give some successful ones, and the inference results for p is correct.

It's actually no related to the bug fixed in that PR but the it's mainly because the choice of priors. Truncated distributions are very sensitive to intializations. In Turing.jl the intializations is done by draw a number from Uniform(-e,e) and transform it to the truanted range (see https://github.com/yebai/Turing.jl/blob/master/src/transform.jl#L35), which works fine for most of the case.

I don't think there is a universal good intialization mechanism for truncated distributions. If they are commonly used in this package, I guess it's better for us to provide an interface for customized initializations for the priors. I can do this and make a PR if it is something wanted.

[ChrisRackauckas]

I don't think there is a universal good intialization mechanism for truncated distributions. If they are commonly used in this package, I guess it's better for us to provide an interface for customized initializations for the priors. I can do this and make a PR if it is something wanted.

That would be very useful. DynamicHMC.jl has a way to specify continuous domain transformations for this purpose. It would be helpful if Turing could read the domain of the prior and directly do a good transformation.

[Vaibhavdixit02]

@xukai I also think it would be a great addition to Turing if it could be done and would be very useful in our case. Also please inform me if I can lend any support, I'll be very glad to be of any help.

[xukai92]

@ChrisRackauckas

That would be very useful. DynamicHMC.jl has a way to specify continuous domain transformations for this purpose. It would be helpful if Turing could read the domain of the prior and directly do a good transformation.

Turing.jl has such a transformation process for variables with constraints, e.g. for Truncated(Normal(10,2),0,15) we always make initlaize it between 0 and 15 according the domain from the prior.

I think here the thing is that when the variables following the truncated Normal are initialized in some region, there are some numerical issue of AD when differentiating some
related functions, which then gives NaN.

I further investigate where the NaN comes from. (Note that when performing HMC we need to do it on unconstrained space, i.e R. For this reason, a variable 'v' in [a, b] is transformed to R by logit((x - a) ./ (b - a)), changed by HMC and transformed back to [a,b] by (b - a) .* invlogit(x) + a), where invlogit{T<:Real}(x::Union{T,Vector{T},Matrix{T}}) = one(T) ./ (one(T) + exp.(-x)) and logit{T<:Real}(x::Union{T,Vector{T},Matrix{T}}) = log(x ./ (one(T) - x)).

1. For some initializations, the initial gradient is very large so the variable in R is changed to something like -1000.
2. When it is transformed back to [a,b], there is numerical problem in AD.

A minimum case for the numerical problem is:

using ForwardDiff: Dual
@inline invlogit{T<:Real}(x::T) = one(T) ./ (one(T) + exp(-x))

d = Dual(-1000, 1)
invlogit(d) # => Dual{Void}(0.0,NaN)

Any idea to improve the stability of the related functions?

[Vaibhavdixit02]

When it is transformed back to [a,b], there is numerical problem in AD.

Can't the AD step be done before the transformation and transformation applied at some point later in the algorithm? I am not very familiar with the details of how HMC is implemented yet so I apologize if this is a bit obtuse on my part.

[xukai92]

I don't think it's possible to do that.

We plan to improve this stability of Turing.jl - let me put this issue in mind when resolving related issues in Turing.jl (TuringLang/Turing.jl#324).

[ChrisRackauckas]

Fixed by #48