I am trying to infer the parameters (a11, a12, …, c21, c22) of the following model:
@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:
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:
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:
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:
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
d = [c for c in eachrow(y)]
and
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?