Turing with multiple replicates of data

I have a simple model for which I have three replicates (data) to be used in inference. I’m new to julia. Right now I’m setting up inference using turing.jl .

My Turing model looks like this :

@model function fit_(data1, data2, data3, prob)
# Prior distributions.
α ~ Uniform(0,0.1)
β ~ Uniform(0,0.1)
γ ~ Uniform(0,0.1)
δ ~ Uniform(0,0.1)
λ ~ Uniform(0,0.2)
ζ ~ Uniform(0,0.2)
σ ~ Uniform(0,2000)

p = [ α, β, γ, δ, λ, ζ]

prob = remake(problem, p=p)

predicted = solve(prob, Tsit5(), saveat=1)'



for i in 1:6

  data1[:,i] .~ MvNormal(predicted[:,i], σ^2 * I)
  data2[:, i] .~ MvNormal(predicted[:,i], σ^2 * I)
  data3[:, i] .~ MvNormal(predicted[:,i], σ^2 * I)
    
end  

return nothing
end

model = fit_(z1, z2, z3 , problem)

Do you guys think its a correct approach ? because model is simple and I inferred same model in r but here i’m not getting correct fitting except for first state in my model.

Model:

function new_model(du, u, p, t)
x, y, z, m, n, r, s = u
α, β, γ, δ, λ, ζ = p
du[1] = - (α + β + γ + δ + λ + ζ ) * x * s
du[2] = α * x * s
du[3] = β * x * s
du[4] = γ * x * s
du[5] = δ * x * s
du[6] = λ * x * s
du[7] = 0
end

true parameters these are random parameters for setting up problem

p = [0.5, 0.2, 0.1, 0.3, 0.4, 0.5]

initial state values

u0 = [2000, 0, 0, 0, 0, 0, 5]
problem = ODEProblem(new_model, u0, (0.0,4.0), p)