# Bayesian regression with parametrized basis functions in Turing.jl

**URL:** https://discourse.julialang.org/t/bayesian-regression-with-parametrized-basis-functions-in-turing-jl/71725
**Category:** Probabilistic Programming
**Tags:** question, turing, gaussian-process
**Created:** [November 18, 2021, 1:49pm UTC](https://discourse.julialang.org/t/bayesian-regression-with-parametrized-basis-functions-in-turing-jl/71725 "2021-11-18T13:49:17Z")
**Posts on this page:** 8
**Page:** 1

<div class="post-metadata">

### Author: ![mvsoom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mvsoom/32/31051_2.png) [@mvsoom](https://discourse.julialang.org/u/mvsoom)
#### Post date: [November 18, 2021, 1:49pm UTC](https://discourse.julialang.org/t/bayesian-regression-with-parametrized-basis-functions-in-turing-jl/71725/1 "2021-11-18T13:49:17Z")

</div>

I’m trying to migrate a simple Bayesian model handwritten in Python to Turing.jl for a more modern approach.

The model is a classical parametrized basis function expansion: a 1D (speech) signal f(t) is expanded into parametrized basis functions, i.e. f(t) = \sum\_k w\_k \phi\_k(t; \theta), where w\_k are the expansion coefficients and \theta are the model parameters. The number of parameters \theta of these basis functions vary from dozens to hundreds, depending on the signal’s complexity and computing budget.

The weights w\_k have a multivariate normal prior and f(t) is assumed to be corrupted by white Gaussian noise. Therefore, the weights can be analytically integrated out of the posterior for \theta and the only parameters that need to be inferred are the \theta.

The catch is that the basis functions \phi\_k(t; \theta) are not known in closed form: they depend on the \theta in a rather complicated way. (They arise from the Karhunen-Loève expansion of a Gaussian process kernel, which depends on the diagonalization of a 2D matrix of Fourier coefficients obtained from the kernel in my case). I am not sure if automatic differentiation would work in this case.

My aim is to perform inference both with optimization (MAP, annealing) and sampling (nested sampling, HMC) approaches. Given my description, do you think that Turing.jl is the appropriate tool for the job?

---

<div class="post-metadata">

### Author: ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [November 18, 2021, 2:00pm UTC](https://discourse.julialang.org/t/bayesian-regression-with-parametrized-basis-functions-in-turing-jl/71725/2 "2021-11-18T14:00:40Z")

</div>

It sounds like you will have to perform a K-L expansion at every step. This is doable but might be more computational complexity than just using a fixed basis set and letting the coefficients be parameters. I think Turing should handle this for you but you might need to try a variety of models to see what is actually efficient.

---

<div class="post-metadata">

### Author: ![mvsoom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mvsoom/32/31051_2.png) [@mvsoom](https://discourse.julialang.org/u/mvsoom)
#### Post date: [November 18, 2021, 5:08pm UTC](https://discourse.julialang.org/t/bayesian-regression-with-parametrized-basis-functions-in-turing-jl/71725/3 "2021-11-18T17:08:02Z")

</div>

The \theta parameters are kind of essential (they have physical interpretations), so I unfortunately cannot use fixed basis functions. The expansion order of K-L is small compared to the number of data points, so it’s not that bad 🙂.

Could you confirm that Turing will recognize automatically that the \theta are parameters to be inferred, even when their likelihood is a complicated nondifferentiable function, given that the `@model` will contain the necessary statements such as `θ₁ ~ Normal(0,1)`?

---

<div class="post-metadata">

### Author: ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [November 18, 2021, 6:41pm UTC](https://discourse.julialang.org/t/bayesian-regression-with-parametrized-basis-functions-in-turing-jl/71725/4 "2021-11-18T18:41:38Z")

</div>

Approximately, things on the left hand side of a Turing ~ statement are parameters unless they appear in the argument to the model function.

```julia
@model function foo(a,b,c)
   theta ~ Normal(0.0,1.0) # theta is a parameter because doesn't appear in arguments to foo
   a ~ Beta(thingy(theta),thingy(b)) # a is data because it's in the arguments to foo, this is a likelihood term
   c ~ CustomDistribution(theta,a,b) # more likelihood
end

```

So Turing should do the right thing.

---

<div class="post-metadata">

### Author: ![ElOceanografo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eloceanografo/32/624_2.png) [@ElOceanografo](https://discourse.julialang.org/u/ElOceanografo)
#### Post date: [November 18, 2021, 7:17pm UTC](https://discourse.julialang.org/t/bayesian-regression-with-parametrized-basis-functions-in-turing-jl/71725/5 "2021-11-18T19:17:15Z")

</div>

I am currently using Turing to model acoustic scattering in the ocean, which involves calculating parameterized basis functions (in my case, frequency-dependent acoustic target strengths). I don’t have as many parameters as you do, but it’s a structurally similar problem and its working for me. Simplified hypothetical example:

```julia
@model function BackscatterExample(freqs, backscatter)
    n ~ Exponential(10) # prior for numerical density of bubbles
    bubble_radius ~ Uniform(0.001, 0.01) # prior for bubble radius
    noise ~ Exponential(1.0) # observation noise
    backscattering_cross_section = map(f -> complicated_scattering_function(bubble_radius, f), freqs)
    pred_backscatter = n * backscattering_cross_section
    backscatter ~ MvNormal(pred_backscatter, noise)
end

```

With a more complicated model performance may become an issue, but at least in principle you should be able to get it to work.

---

<div class="post-metadata">

### Author: ![mvsoom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mvsoom/32/31051_2.png) [@mvsoom](https://discourse.julialang.org/u/mvsoom)
#### Post date: [November 19, 2021, 9:20am UTC](https://discourse.julialang.org/t/bayesian-regression-with-parametrized-basis-functions-in-turing-jl/71725/6 "2021-11-19T09:20:22Z")

</div>

Thank you both very much for your answers. I will start implementing today.

Cheers!

---

<div class="post-metadata">

### Author: ![theogf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/theogf/32/1987_2.png) [@theogf](https://discourse.julialang.org/u/theogf)
#### Post date: [April 19, 2022, 11:39am UTC](https://discourse.julialang.org/t/bayesian-regression-with-parametrized-basis-functions-in-turing-jl/71725/7 "2022-04-19T11:39:07Z")

</div>

A bit late to the party but you might be interested in [https://github.com/JuliaGaussianProcesses/BayesianLinearRegressors.jl](https://github.com/JuliaGaussianProcesses/BayesianLinearRegressors.jl)

---

<div class="post-metadata">

### Author: ![mvsoom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mvsoom/32/31051_2.png) [@mvsoom](https://discourse.julialang.org/u/mvsoom)
#### Post date: [April 19, 2022, 5:55pm UTC](https://discourse.julialang.org/t/bayesian-regression-with-parametrized-basis-functions-in-turing-jl/71725/8 "2022-04-19T17:55:14Z")

</div>

Yes, I am using that library :). Thanks for the suggestion!
