Hi @svilupp!
The way you generate data is the following, right?
# generate data
a=0.2
b=5
time_index=collect(1.:20.)
noise_scale=0.5
y=time_index .* a .+ b + noise_scale*randn(size(time_index,1))
Here, the noise_scale variable corresponds to the standard deviation 0.5 (or variance 0.25).
Note that the sigma variable of your model corresponds to the precision parameter, not the variance.
As you might know, the relationship between these variables is the following: 1/variance = precision.
When you retrieve the posterior of the precision parameter and check its mean and variance (corresponding to the precision of your likelihood) you get:
msigma = mean(results.posteriors[:sigma]) # mean of the posterior
vsigma = var(results.posteriors[:sigma]) # variance of the posterior
to retrieve the estimates for noise_scale, we need to apply inversion and take a square root of precision estimates:
inv_msigma = sqrt(inv(msigma)) # estimates for noise_scale
In this way, I get inv_msigma = 0.49521700118560297!
You can try to generate the data in the following way (so it maps to your model likelihood ±one-to-one):
y = time_index .* a .+ b .+ rand(NormalMeanPrecision(0.0, noise_precision), size(time_index,1));