# BoundsError when inverting a model with Turing

**URL:** <https://discourse.julialang.org/t/boundserror-when-inverting-a-model-with-turing/74019>\
**Category:** New to Julia\
**Tags:** question, diffeq, turing, arrays\
**Created:** [January 4, 2022, 8:46am UTC](https://discourse.julialang.org/t/boundserror-when-inverting-a-model-with-turing/74019 "2022-01-04T08:46:42Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![spinystellate](https://avatars.discourse-cdn.com/v4/letter/s/e8c25b/32.png) [@spinystellate](https://discourse.julialang.org/u/spinystellate)\
**Post date:** [January 4, 2022, 8:46am UTC](https://discourse.julialang.org/t/boundserror-when-inverting-a-model-with-turing/74019/1 "2022-01-04T08:46:42Z")

</div>

I am trying to infer the parameters (a11, a12, …, c21, c22) of the following model:

```julia
@model function invert(u, y, pE, pC, options)
    
    # Specify priors
    h = [Normal(pE.H[i], pC.H[i]) for i in eachindex(pE.H)]
    hc = [Normal(pE.HC[i], pC.HC[i]) for i in eachindex(pE.HC)]
    H = reshape(rand.(h), 2, 7)
    HC = reshape(rand.(hc), 2, 2)
    
    a11 ~ Normal(pE.A[1,1], pC.A[1,1])
    a12 ~ Normal(pE.A[1,2], pC.A[1,2])
    a21 ~ Normal(pE.A[2,1], pC.A[2,1])
    a22 ~ Normal(pE.A[2,2], pC.A[2,2])
    
    c11 ~ Normal(pE.C[1,1], pC.C[1,1])
    c12 ~ Normal(pE.C[1,2], pC.C[1,2])
    c21 ~ Normal(pE.C[2,1], pC.C[2,1])
    c22 ~ Normal(pE.C[2,2], pC.C[2,2])
    
    x0 = zeros(options.ℓ, 1) # initial condition
    z0 = zeros(options.ℓ, 6) # initial condition
    states = [x0, z0]
    P = [a11, a12, a21, a22, c11, c12, c21, c22, u, H, HC]
    tspan = (0.0, options.duration)
    predicted = forward_model(P, states, options)  
 
    y[:,1] ~ Normal(predicted[:,1]) 
    y[:,2] ~ Normal(predicted[:,2])
 
end

m = invert(U, y, pE, pC, options)
chain = mapreduce(c -> sample(m, NUTS(.65),1000), chainscat, 1:3) 

```

but get the error message:

`BoundsError: attempt to access 1-element Vector{Matrix{Float64}} at index [2]`

The error occurs in a function called g, which is called by the function forward\_model:

```julia
function forward_model(params, states, options)
    a11, a12, a21, a22, c11, c12, c21, c22, u, H, HC = params
    x0, z0 = states
    
    ts = range(0.05, step=options.dt, stop=options.duration) # sample times
    tspan = (0.0, options.duration) # time span
    
    p = [a11, a12, a21, a22, c11, c12, c21, c22, u, options.dt] # parameters
    ode = ODEProblem(f, x0, tspan, p)
    sol = solve(ode, Tsit5(), saveat=ts, tstops=ts, dt=1e-2)
    
    ph = [H, HC, sol.u, options] # parameters
    odeh = ODEProblem(h, z0, tspan, ph)
    solh = solve(odeh, Tsit5(), saveat=ts, tstops=ts, dt = 1e-2)
    
    # Model signal
    y = g(solh.u, H, options)
end

function g(z, H, opt)
    # Define constants
    TE = 0.04
    V₀ = 100*0.04
    r₀ = 110
    ϑ₀ = 84.8
    E₀ = opt.Pf_H[:,5] .* exp.(H[:,5]) 
    εₕ = opt.Pf_H[2,6] .* exp.(H[2,6])
    
    # Compute coefficients
    k₁ = 4.3.*ϑ₀.*E₀.*TE
    k₂ = εₕ.*r₀.*E₀.*TE
    k₃ = 1 - εₕ
    
    # Output equation
    v = map(i -> (z[i][:,3]), 1:12000)
    v = transpose(reduce(hcat, v))
    q = map(i -> (z[i][:,4]), 1:12000)
    q = transpose(reduce(hcat, q))
    y = V₀ * (k₁'.*(1 .- q) + k₂'.*(1 .- q./v) + k₃'.*(1 .- v))
end

```

For the sake of simplicity I have left out all other functions (f, h, etc.); I think that these should not be necessary for answering my question. The input z into function g is the solution of an ODE problem and has the following form:

```julia
12000-element Vector{Matrix{Float64}}:
 [0.0 1.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0]
 [0.0 1.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0]
 ⋮

```

Generating data with the function forward\_model works. The BoundsError during the inference procedure occurs in the line:

```julia
v = map(i -> (z[i][:,3]), 1:12000)

```

(in function g). I can see that the error occurs because the input argument y (in function invert) is a Matrix{Float64} of size (12000, 2), while z (in function g) is a Vector{Matrix{Float64}} of size (12000,), but my attempts to resolve the problem remain unsuccessful.

I have tried, for example, to convert the data in the function invert:

```julia
d = [c[:] for c in eachrow(y)]
d[:,1] ~ Normal(predicted[:,1]) 
d[:,2] ~ Normal(predicted[:,2])

```

but this converts d to a Vector{Vector{Float64}} and leads to the same BoundsError. I also tried

```julia
d = [c for c in eachrow(y)]

```

and

```julia
d = [(y)]

```

but these are also not correct. How can I convert y to a Vector{Matrix{Float64}}? Or is there a better alternative solution?

---

<div class="post-metadata">

**Author:** ![spinystellate](https://avatars.discourse-cdn.com/v4/letter/s/e8c25b/32.png) [@spinystellate](https://discourse.julialang.org/u/spinystellate)\
**Post date:** [January 4, 2022, 3:03pm UTC](https://discourse.julialang.org/t/boundserror-when-inverting-a-model-with-turing/74019/2 "2022-01-04T15:03:58Z")

</div>

Ok, the following converts the data to a Vector{Matrix{Float64}} of size (12000,):

```julia
d = [hcat(c[:]) for c in eachrow(data)]
d[:,1] ~ Normal(predicted[:,1]) 
d[:,2] ~ Normal(predicted[:,2])

```

but this was not the solution after all - I still get the same BoundsError in the same place.

Any suggestions of what else may be wrong with my code/approach would be very welcome.