# Parameterizing matrix in Turing model without mutating arrays

**URL:** <https://discourse.julialang.org/t/parameterizing-matrix-in-turing-model-without-mutating-arrays/56977>\
**Category:** New to Julia\
**Tags:** turing\
**Created:** [March 11, 2021, 9:49pm UTC](https://discourse.julialang.org/t/parameterizing-matrix-in-turing-model-without-mutating-arrays/56977 "2021-03-11T21:49:16Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![parb](https://avatars.discourse-cdn.com/v4/letter/p/4da419/32.png) [@parb](https://discourse.julialang.org/u/parb)\
**Post date:** [March 11, 2021, 9:49pm UTC](https://discourse.julialang.org/t/parameterizing-matrix-in-turing-model-without-mutating-arrays/56977/1 "2021-03-11T21:49:16Z")

</div>

Hi All.

So I’ve got this compartment model which uses a matrix exponential,  
\hat y = e^{A\,t} x\_0,  
where the elements of the N by N matrix A are flow rates between N compartments, and x\_0 are the initial compartment concentrations. I hope to estimate these flow rates in A.

The matrix A should be only positive numbers on the off-diagonals, and the diagonals should be their negative column sums. How do I achieve this without mutating arrays?

I can produce A using the Zygote backend by producing a multivariate normal and then reshaping it, as I’ve done below (for the case of 4 compartments).

Is there a way to parameterize matrices without mutating arrays?

Or should I perhaps implement the ODE behind this model instead?

Thanks a bunch for checking my question.

Here’s what I _can_ do:

```julia
@model function makemodel(data, init, time_range)
    N = size(data)[1]
    σ1 ~ InverseGamma(2, 2)
    σ2 ~ InverseGamma(10, 9)
    σ3 ~ InverseGamma(3, 0.5)
    Sigma = zeros(16, 16) + Diagonal([σ2, σ1, σ1, σ1,
                                      σ1, σ2, σ1, σ1,
                                      σ1, σ1, σ2, σ1,
                                      σ1, σ1, σ1, σ2])
    K ~ MvNormal(Sigma)
    A = reshape(K, N, N)
    
    predicted = reduce(hcat, [real(exp(A*t))*init for t in time_range])
    
    data = transpose(data)[:]
    data ~ arraydist([Normal(predicted[i], σ3) for i=1:length(predicted)])
end

```

Here’s what I _wish_ I could do, for the case of 3 compartments:

```julia
@model function makemodel(data, init, time_range)
    N = 3
    σ1 ~ InverseGamma(2, 2)
    σ2 ~ InverseGamma(3, 0.5)
    K ~ MvNormal(3, σ2)
    K_diags ~ TruncatedNormal([0.0, 0.0, 0.0], 1.0, -100.0, 0.0)
    A = zeros(3, 3)
    A[2, 1] = K[1]
    A[3, 1] = K[2]
    A[3, 2] = K[2]
    for i=1:3
        A[i, i] = -sum(A[:, i]) + K_diags[i]
    
    predicted = reduce(hcat, [real(exp(A*t))*init for t in time_range])
    
    data = transpose(data)[:]
    data ~ arraydist([Normal(predicted[i], σ3) for i=1:length(predicted)])
end

```