Hi folks,

I was thinking about this implementation to sample `$\epsilon$`

with multivariate distribution. Similar to the way we handle variables in random block. Here is an example: (See comments below the code):

```
const N_THETA = 3
const N_OMEGA = 3
const N_SIGMA = 2
model = @model begin
@param begin
θ ∈ VectorDomain(N_THETA)
Ω ∈ PDiagDomain(N_OMEGA)
Σ ∈ PDiagDomain(N_SIGMA)
end
@random begin
η ~ MvNormal(Ω)
ϵ ~ MvNormal(Σ)
end
@pre begin
CLᵢ = exp(log(θ[1]) + η[1])
Vᵢ = exp(log(θ[2]) + η[2])
kaᵢ = exp(log(θ[3]) + η[3])
keᵢ = CLᵢ/Vᵢ
end
@dynamics begin
y1' = - kaᵢ * y1
y2' = kaᵢ* y1 - keᵢ * y2
end
@derived begin
cp = @. y2/Vᵢ
dv = @. cp + cp * ϵ[1] + ϵ[2]
end
end
```

**COMMENT**: Currently, epsilon will have one value per subject and will not change which is expected from `@random`

block. Is there a way that one can sample within the same subject?

Here is a dummy setup:

```
event = DosageRegimen(100)
sub = Subject(;id=1, events=event)
par_values = (θ=[8, 25, 1.5], Ω = Diagonal([0.09, 0.09, 0.09]), Σ = Diagonal([0.015, 0.0625]), )
sim = simobs(model, sub, par_values) |> DataFrame
```

You can see the values of epsilons from dataframe.