# Using DynamicHMC for a model defined in Python

**URL:** <https://discourse.julialang.org/t/using-dynamichmc-for-a-model-defined-in-python/117765>\
**Category:** Specific Domains\
**Tags:** dynamichmc\
**Created:** [August 2, 2024, 12:47pm UTC](https://discourse.julialang.org/t/using-dynamichmc-for-a-model-defined-in-python/117765 "2024-08-02T12:47:15Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![Douwe](https://avatars.discourse-cdn.com/v4/letter/d/7ea924/32.png) [@Douwe](https://discourse.julialang.org/u/Douwe)\
**Post date:** [August 2, 2024, 12:47pm UTC](https://discourse.julialang.org/t/using-dynamichmc-for-a-model-defined-in-python/117765/1 "2024-08-02T12:47:15Z")

</div>

Hi I’ve been trying to use **DynamicHMC** to sample my model written in **Python**. The problem I’m encountering is that I don’t know how to structure the output of my **log posterior** function in julia correctly, so the MCMC sampler knows what to do with it.

In the examples provided it is suggested to use the following code  
t = problem\_transformation(p)  
P = TransformedLogDensity(t, p)  
∇P = ADgradient(:ForwardDiff, P);  
[https://www.tamaspapp.eu/DynamicHMCExamples.jl/latest/example\_linear\_regression/](https://www.tamaspapp.eu/DynamicHMCExamples.jl/latest/example_linear_regression/)

However, because I use python output, the **gradient** of the likelihood is not calculated automatically anymore. So currently I only compute the log posterior without the gradients. And even if I would calculate these gradients myself. I am not sure how to give this all the proper **structure** so the sampler mcmc\_with\_warmup() can use it.

For context the function that computes the logposterior in the example linked above returns output such as (I printed it):  
Dual{ForwardDiff.Tag{Base.Fix1{typeof(LogDensityProblems.logdensity), TransformedLogDensity{TransformVariables.TransformTuple{@NamedTuple{β::TransformVariables.ArrayTransformation{TransformVariables.Identity, 1}, σ::TransformVariables.ShiftedExp{true, Int64}}}, LinearRegressionProblem{Vector{Float64}, Matrix{Float64}, Float64}}}, Float64}}(-82.04794174652842,27.668542512145063,-19.965222890254054,-5.894398700546053,10.918972471706919)

This is my julia code:

```julia
ENV["JULIA_CONDAPKG_BACKEND"] = "Null" 
import PythonCall
py_module = PythonCall.Core.pyimport("models_douwe") 
using TransformVariables, TransformedLogDensities, LogDensityProblems, LogDensityProblemsAD,
    DynamicHMC, DynamicHMC.Diagnostics, SimpleUnPack, Statistics, Random, Distributions, 
    ThreadTools, MCMCDiagnosticTools, BenchmarkTools

struct Model_Structure
    y::Vector # hydraulic head observations
    model_nr::Int # identifier of the model defined in FloPy
    locations::Tuple # measurement locations
end

function create_model_structure(model_nr::Int, locations::Tuple)
    y_py = py_module.observations(model_nr, locations)
    y_jl = PythonCall.Core.pyconvert(Tuple, y_py)
    y = collect(y_jl)
    p = Model_Structure(y, model_nr, locations)
    return p
end

function (problem::Model_Structure)(θ) 
    @unpack hk, vk = θ 
    @unpack y, model_nr, locations = problem
    #create input suitable for python 
    if typeof(hk) == Float64
        all_k = (hk, vk) 
    else
        all_k = (hk.value, vk.value) 
    end
    #calculate likelihood
    model_output_py = py_module.model_output(model_nr, locations, K = all_k) 
    model_output_jl = PythonCall.pyconvert(Tuple, model_output_py) # julia Tuple
    model_output = collect(model_output_jl) # hydrualic head measurements, as a vector
    ϵ_distribution = Normal(0, 0.1)
    ℓ_error = loglikelihood(ϵ_distribution, y - model_output)
    # calculate prior density
    ℓ_hk = loglikelihood(Normal(10, 10), hk) 
    ℓ_vk = loglikelihood(Normal(10, 10), vk)
    return ℓ_error + ℓ_hk + ℓ_vk
end

function problem_transformation(p::Model_Structure)
    as((hk = asℝ₊, vk = asℝ₊)) # hydraulic conductivity cannot be negative
end

p = create_model_structure(1, LOCATIONS) # redefine for clarity
trans = problem_transformation(p) #transforms all k's to real positive numbers
P = TransformedLogDensity(trans, p) 
∇P = ADgradient(:ForwardDiff, P)

results = mcmc_with_warmup(Random.default_rng(), ∇P, 1000)

```

Any suggestions on correctly incorporating the gradients in my case?