# Gaussian process with noisy heteroscedastic correlated targets (outputs)

**URL:** <https://discourse.julialang.org/t/gaussian-process-with-noisy-heteroscedastic-correlated-targets-outputs/132582>\
**Category:** Statistics\
**Created:** [September 22, 2025, 7:31pm UTC](https://discourse.julialang.org/t/gaussian-process-with-noisy-heteroscedastic-correlated-targets-outputs/132582 "2025-09-22T19:31:06Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![mocalvao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mocalvao/32/19318_2.png) [@mocalvao](https://discourse.julialang.org/u/mocalvao)\
**Post date:** [September 22, 2025, 7:31pm UTC](https://discourse.julialang.org/t/gaussian-process-with-noisy-heteroscedastic-correlated-targets-outputs/132582/1 "2025-09-22T19:31:06Z")

</div>

Hi there,  
I have 2 kinds of datasets to which I would like to fit a Gaussian process:  
(1) a dataset with N observations: (x\_i, y\_i, \sigma\_i) , i=1, 2, ..., N, where x\_i are the input (covariate) variables, y\_i the output (target) variables, and \sigma\_i are the varying noises, so that it is a heteroscedastic process. I was able to use, in Python, the `sci-kit learn`’s package GaussianProcessRegressor, with the parameter `alpha` equal to the vector \sigma and everything went ok. I guess there should be an equally easy procedure in one of the packages of the GaussianProcesses echosystem of Julia. Could you give me the exact reference?  
(2) another dataset, constituted by two files:

- the first one is the observations (x\_i, y\_i) , i=1, 2,...N
- the second one contains the noise for the output variables (targets) in the form of a N \times N covariance matrix (array) \Sigma , no longer diagonal. This is what I have called noisy heteroscedastic correlated targets in the title of my topic. Does it make sense to fit a Gaussian process to such a dataset? If so, which Julia package allows implementing this and exactly how? Pondering a bit deeper, I think (is this correct?) what this second file amounts to is a sample or estimate for the `covariance function` k(x\_i, x\_j) , i,j =1, 2, ..., N ; as such I should use some approximation to it to define a new kernel for my call to the Gaussian process routine…

Thank you very much in advance!

---

<div class="post-metadata">

**Author:** ![eteppo](https://avatars.discourse-cdn.com/v4/letter/e/90db22/32.png) [@eteppo](https://discourse.julialang.org/u/eteppo)\
**Post date:** [September 25, 2025, 8:52am UTC](https://discourse.julialang.org/t/gaussian-process-with-noisy-heteroscedastic-correlated-targets-outputs/132582/2 "2025-09-25T08:52:11Z")

</div>

I can point out one possible approach, using Turing.jl. For example, your model could look like below?

```julia-auto
using Turing, LinearAlgebra
@model function GPR(x, Σ)
    N = length(x)
    η ~ Exponential(2)
    ρ ~ Exponential(1)
    d² = (x .- x').^2
    Κ = η * (exp.(-ρ * d²))
    y ~ MvNormal(zeros(N), Κ + Σ)
end

```

No guarantees about anything (don’t know anything about GPs), just to note that it lets you try almost anything. ([Turing Gaussian Process tutorial](https://turinglang.org/docs/tutorials/gaussian-processes-introduction/))

---

<div class="post-metadata">

**Author:** ![mocalvao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mocalvao/32/19318_2.png) [@mocalvao](https://discourse.julialang.org/u/mocalvao)\
**Post date:** [September 25, 2025, 11:12am UTC](https://discourse.julialang.org/t/gaussian-process-with-noisy-heteroscedastic-correlated-targets-outputs/132582/3 "2025-09-25T11:12:43Z")

</div>

Thanks @eteppo,  
I will ponder about your suggestion.