# HPV disease model estimation (beginner issue)

`Hi, I 'm trying to estimate HPV model parameters with Turing but I’m not sure about how to determine prior parameters especially with σ , so when I run the model it seems because of the bad selection of σ last parameter takes high values (it has to be between 0-1 but it seems chain stars near 1000.0 values for the last parameter).Lastly I want to compose my observation vector with two prevalence values(one for females and one for males) and two incidence values , so I guess I had to determine two different σ values, but I m not sure about σ has to be univariate or multivariate distribution.
Sorry for beginner code but I do not have advisor in my university , this will be my first attempt for feedback. I’m really stuck with that issue and I don’t even know where I lack knowledge any help will make my day.
Thanks.
Mevlüt.

``````using DifferentialEquations, LinearAlgebra, ModelingToolkit, Plots, Parameters, Turing, Distributions
using MCMCChains, StatsPlots
using Random

Random.seed!(14000);

using Logging
Logging.disable_logging(Logging.Warn)

const J = 3
const K = 5
const TF = 20
const sriskgr = range(0,stop=1,length=J)
const agegr = range(0,stop=1,length=K)
const timefinale = range(0,stop=1,length=TF)
const NF = 1500000
const NM = 1500000
const Rs = [6/10,3/10,1/10]
const Ra = [3/10,1/4,1/5,3/20,1/10]
const prv = 0.043

function init_nfemale(js,ka)
J = length(js)
K = length(ka)
Nfemale = zeros(J,K)
rsa = Rs*transpose(Ra)
Nfemale = NF*rsa
Nfemale
end

const Nfemale = init_nfemale(sriskgr,agegr)

function init_nmale(js,ka)
J = length(js)
K = length(ka)
Nmale = zeros(J,K)
rsa = Rs*transpose(Ra)
Nmale = NM*rsa
Nmale
end

const Nmale = init_nmale(sriskgr,agegr)

function init_smat(js,ka)
J = length(js)
L = length(js)
K = length(ka)
M = length(ka)
cstar = 2

actc = zeros(J,L,K,M)
actcmult = zeros(J,K)
actc[1,1,:,:] = I(K)
actc[2,2,:,:] = I(K)
actc[3,3,:,:] = I(K)
actcnew = mapslices(sum, actc, dims = [2,4])
actcnew2 = dropdims(actcnew; dims=2)
actcnew3 = dropdims(actcnew2; dims=3)
actcmult =  cstar * actcnew3

actcmult
end

const actcmult = init_smat(sriskgr,agegr)

function ODEHPV(du,u,p,t)
DIMM, DIP, DCN, DCC, DIPm, PCI, PCCN, PCC, BTRM, BTRF  = p

@inbounds for j = 1, i in 1:J
du[i,j,1] = u[i,j,5]/DIMM - u[i,j,1]*(BTRF*actcmult[i,j]*u[i,j,7]/Nmale[i,j]) - 0.1*u[i,j,1] + 0.1*Rs[i]*(u[5,j,2]+u[5,j,3]+u[5,j,4]+u[5,j,5])
du[i,j,2] = u[i,j,1]*(BTRF*actcmult[i,j]*u[i,j,7]/Nmale[i,j]) - u[i,j,2]/DIP - 0.1*u[i,j,2]
du[i,j,3] = (1-PCI)*u[i,j,2]/DIP - u[i,j,3]/DCN - 0.1*u[i,j,3]
du[i,j,4] = (1-PCCN)*u[i,j,3]/DCN - u[i,j,4]/DCC - 0.1*u[i,j,4]
du[i,j,5] = PCI*u[i,j,2]/DIP + PCCN*u[i,j,3]/DCN + PCC*u[i,j,4]/DCC - u[i,j,5]/DIMM - 0.1*u[i,j,5]
du[i,j,6] = u[i,j,8]/DIMM - u[i,j,6]*(BTRM*actcmult[i,j]*u[i,j,2]/Nfemale[i,j]) - 0.1*u[i,j,6] + 0.1*Rs[i]*(u[5,j,7]+u[5,j,8])
du[i,j,7] = u[i,j,6]*(BTRM*actcmult[i,j]*u[i,j,2]/Nfemale[i,j]) - u[i,j,7]/DIPm - 0.1*u[i,j,7]
du[i,j,8] = u[i,j,7]/DIPm - u[i,j,8]/DIMM - 0.1*u[i,j,8]

end

@inbounds for j in 2:K, i in 1:J
du[i,j,1] = u[i,j,5]/DIMM - u[i,j,1]*(BTRF*actcmult[i,j]*u[i,j,7]/Nmale[i,j]) - 0.1*u[i,j,1] + 0.1*u[i,j-1,1]
du[i,j,2] = u[i,j,1]*(BTRF*actcmult[i,j]*u[i,j,7]/Nmale[i,j]) - u[i,j,2]/DIP - 0.1*u[i,j,2] + 0.1*u[i,j-1,2]
du[i,j,3] = (1-PCI)*u[i,j,2]/DIP - u[i,j,3]/DCN - 0.1*u[i,j,3] + 0.1*u[i,j-1,3]
du[i,j,4] = (1-PCCN)*u[i,j,3]/DCN - u[i,j,4]/DCC - 0.1*u[i,j,4] + 0.1*u[i,j-1,4]
du[i,j,5] = PCI*u[i,j,2]/DIP + PCCN*u[i,j,3]/DCN + PCC*u[i,j,4]/DCC - u[i,j,5]/DIMM - 0.1*u[i,j,5] + 0.1*u[i,j-1,5]
du[i,j,6] = u[i,j,8]/DIMM - u[i,j,6]*(BTRM*actcmult[i,j]*u[i,j,2]/Nfemale[i,j]) - 0.1*u[i,j,6] + 0.1*u[i,j-1,6]
du[i,j,7] = u[i,j,6]*(BTRM*actcmult[i,j]*u[i,j,2]/Nfemale[i,j]) - u[i,j,7]/DIPm - 0.1*u[i,j,7] + 0.1*u[i,j-1,7]
du[i,j,8] = u[i,j,7]/DIPm - u[i,j,8]/DIMM - 0.1*u[i,j,8] + 0.1*u[i,j-1,8]
end

end

function init_mevotoar(js,ka)
J = length(js)
K = length(ka)
u = zeros(J, K, 8)
Rs = [6/10,3/10,1/10]
Ra = [3/10,1/4,1/5,3/20,1/10]
rsa = Rs*transpose(Ra)
u[:,:,1] = (1-prv)* (NF * rsa)
u[:,:,2] = prv* (NF * rsa)
u[:,:,3] = zeros(3,5)
u[:,:,4] = zeros(3,5)
u[:,:,5] = zeros(3,5)
u[:,:,6] = (1-prv)* (NM * rsa)
u[:,:,7] = prv* (NM * rsa)
u[:,:,8] = zeros(3,5)
u
end

p = [8.0,6.0,8.0,5.0,5.0,0.88,0.65,0.6,0.1,0.1]

u0 = init_mevotoar(sriskgr,agegr)

tspan = (1.0,20.0)

prob = ODEProblem(ODEHPV,u0,tspan,p)

odep = modelingtoolkitize(prob)

jac = eval(ModelingToolkit.generate_jacobian(odep)[2])

f = ODEFunction(ODEHPV, jac=jac)

prob_jac = ODEProblem(f,u0,tspan,p)

serodataf = fill(0.043,(1,TF))
serodatam = fill(0.043,(1,TF))
datacnstage = fill(1430.0,(1,TF))
dataccstage = fill(85.0,(1,TF))

@model function fitlv2(data,prob_jac) # data should be a Vector
σ ~ InverseGamma(2, 3)
DIMM ~ Uniform(4.0,12.0)
DIP ~ Uniform(4.0,8.0)
DCN ~ Uniform(6.0,10.0)
DCC ~ Uniform(4.0,6.0)
DIPm ~ Uniform(2.0,8.0)
PCI ~ Beta(2,2)
PCCN ~ Beta(2,2)
PCC ~ Beta(2,2)
BTRM ~ Beta(2,5)
BTRF ~ Beta(2,5)

p = [DIMM,DIP,DCN,DCC,DIPm,PCI,PCCN,PCC,BTRM,BTRF]
prob = remake(prob_jac, p=p)
predictednonesum = solve(prob,AutoTsit5(Rosenbrock23()), jac = true, saveat = 1.0)
predicted0 = Array(predictednonesum)

sumpredicted = mapslices(sum,predicted0, dims = 1)
sumpredicted1 = dropdims(sumpredicted; dims=1)

predInfectedFemale = Array(sumpredicted1[:,2,:])
predInfectedMale = Array(sumpredicted1[:,7,:])

predCNStagef = Array(sumpredicted1[:,3,:])
predCCStagef = Array(sumpredicted1[:,4,:])

pred_infectedfemale = mapslices(sum,predInfectedFemale, dims = 1)

pred_infectedmale = mapslices(sum,predInfectedMale, dims = 1)

predSeroPrevelancef = similar(pred_infectedfemale)
predSeroPrevelancem = similar(pred_infectedmale)

predSeroPrevelancef = pred_infectedfemale/NF
predSeroPrevelancem = pred_infectedmale/NM

predsumCNStagef = mapslices(sum, predCNStagef, dims = 1)

predsumCCStagef = mapslices(sum, predCCStagef, dims = 1)

predictedfnl = vcat(predSeroPrevelancef,predSeroPrevelancem,predsumCNStagef,predsumCCStagef)

for i = 1:length(predictedfnl)
data[i] ~ Normal(predictedfnl[i], σ) # predicted[i][2] is the data for y - a scalar, so we use Normal instead of MvNormal
end
end

model2 = fitlv2(odedatarndm,prob_jac)

Transform your values so that `[-infty,infty]` maps to `[0,1]`. You can do something like `sigma = 1/(1+exp(-_sigma)`, and then you do parameter estimation on `_sigma` and transform inside of the model. In the end, it’ll estimate the correct `_sigma` and you then just have to transform it to interpret the results.