# Target Variable is a Deterministic Function of Random Variables
**URL:** https://discourse.julialang.org/t/target-variable-is-a-deterministic-function-of-random-variables/71935
**Category:** Probabilistic Programming
**Tags:** turing, distributions, bayesian-inference
**Created:** [November 23, 2021, 10:50am UTC](https://discourse.julialang.org/t/target-variable-is-a-deterministic-function-of-random-variables/71935 "2021-11-23T10:50:26Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@Stephen2](https://discourse.julialang.org/u/Stephen2)
#### Post date: [November 23, 2021, 10:50am UTC](https://discourse.julialang.org/t/target-variable-is-a-deterministic-function-of-random-variables/71935/1 "2021-11-23T10:50:27Z")
Crossposting this from “New to Julia” as I received no responses.
Hi, I am new to Julia and struggling to formulate my model correctly. My model consists of:
- A regressor **x\_t**
- An intermediate variable **y\_t** = **k** \* x\_t + ϵ , ϵ ~ Normal(0, **σ** )
And finally a target variable z that is a function of y\_t and its lags:
- **z\_t** = λ \* y\_t + λ^2 \* y\_(t-1) + λ^3 \* y\_(t-2) + …, **λ** ~ Uniform(0, 1)
The variables of interest are **k, σ** and **λ.**
My issue is that I am unable to tell Julia that z\_t is distributed as above, because by the time I get to defining it, the terms are numbers instead of distributions.
It might be possible to re-parametrize and define z\_t as a normally distributed variable with mean equal to z\_t above, and some variance that is a function of σ and λ. However I’m not willing to attempt to find this function because I’m sure there is a better way.
Below is a version of my code that produces semi-reasonable code, using 4 lags. I’ve attempted to use Distributions.convolute() multiple times to combine distributions, but I’m not sure how correct this is, and again I’m sure there is a better way.
```julia
using Turing, StatsPlots, Random
Random.seed!(123)
function lag_transform(λ, y)
# apply the lag transformation to y
y_type = promote_type(eltype(λ), eltype(y))
z = zeros(y_type, 0)
for i in 1:length(y)
powers = reverse([λ ^ k for k in range(1, stop=i)])
previous_terms = zeros(y_type, 0)
for j in 1:i
append!(previous_terms, y[j] * powers[j])
end
current_z_val = sum(previous_terms)
append!(z, current_z_val)
end
return z
end
# generate some test data
n_samples = 200
σ = 0.5
k = 2.3
λ = 0.7
x = 1:n_samples
y = [x_val * k + rand(Normal(0, σ)) for x_val in x]
z = lag_transform(λ, y)
@model function model_func(x, z)
k ~ Uniform(0.5, 3)
λ ~ Uniform(0.4, 1)
σ ~ truncated(Normal(0, 1), 0, Inf)
y = TArray([convert(Real, k * x_val) for x_val in x])
powers = TArray([convert(Real, λ ^ i) for i in 1:5])
for n in 5:length(y) # start on 5, because we use the previous 4 lags
normals = [
Normal(powers[1] * y[n], powers[1] * σ),
Normal(powers[2] * y[n - 1], powers[2] * σ),
Normal(powers[3] * y[n - 2], powers[3] * σ),
Normal(powers[4] * y[n - 3], powers[4] * σ),
Normal(powers[5] * y[n - 4], powers[5] * σ)
]
z[n] ~ reduce(convolve, normals)
end
end
model = model_func(x, z[5:end]) # dont pass all of z_t, because we use the previous 4 lags
chain = sample(model, NUTS(), 5000)
plot(chain)
```
I’ve also tried some other solutions (looking into @submodels, !addlogprob, etc.) but have been completely unsuccessful. Thanks!
---
### Author:  [@Stephen2](https://discourse.julialang.org/u/Stephen2)
#### Post date: [November 24, 2021, 11:29am UTC](https://discourse.julialang.org/t/target-variable-is-a-deterministic-function-of-random-variables/71935/2 "2021-11-24T11:29:03Z")
Apologies, `z[n] ~ reduce(convolve, normals)` should be `z[n - 4] ~ reduce(convolve, normals)`.
---
### Author:  [@mikkoku](https://discourse.julialang.org/u/mikkoku)
#### Post date: [November 24, 2021, 2:34pm UTC](https://discourse.julialang.org/t/target-variable-is-a-deterministic-function-of-random-variables/71935/3 "2021-11-24T14:34:22Z")
I think the problem is in the model formulation. It looks like z\_t = \lambda z\_{t-1} + \lambda y\_t and this looks like an AR(1) model. If this is the case then it is easy to find the distribution of z\_t, namely z\_t \sim Normal(\lambda (z\_{t-1} + kx\_t), \lambda \sigma).
---
### Author:  [@Stephen2](https://discourse.julialang.org/u/Stephen2)
#### Post date: [November 24, 2021, 4:20pm UTC](https://discourse.julialang.org/t/target-variable-is-a-deterministic-function-of-random-variables/71935/4 "2021-11-24T16:20:14Z")
Thanks for the response. Unfortunately the model isn’t z\_t = \lambda z\_{t-1} + \lambda y\_t. The transformation I want to apply is the [Adstock decay effect](https://en.wikipedia.org/wiki/Advertising_adstock#Advertising_lag:_decay_effect), where the decay function decreases geometrically: z\_t = \lambda y\_t + \lambda ^2y\_{t-1} + \lambda ^ 3 y\_{t-2} + \dots (don’t worry about the lack of normalization, I do this later).
---
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [November 24, 2021, 5:43pm UTC](https://discourse.julialang.org/t/target-variable-is-a-deterministic-function-of-random-variables/71935/5 "2021-11-24T17:43:08Z")
> [@Stephen2](#):
>
> My issue is that I am unable to tell Julia that z\_t is distributed as above, because by the time I get to defining it, the terms are numbers instead of distributions.
Do I understand correctly that you believe `z = f(SomeRandomStuff) + no_error`
You will never get good sampling but another way to put this is
```julia
z ~ Dirac(f(SomeRandomStuff))
```
If this is the case, it’s similar to ABC methods where you have a simulation model and are trying to match data, but you can’t match it exactly so you have to tolerate some error. I’d suggest
```julia
z ~ Normal(f(SomeRandomStuff), epsilon)
```
where epsilon is small, now the size of epsilon is something you can tune to compromise between the Dirac type situation (where you will never get samples at all) and the situation where epsilon is large and you’ll get samples but most of them will be too far off your target to be acceptable.
---
### Author:  [@Stephen2](https://discourse.julialang.org/u/Stephen2)
#### Post date: [November 25, 2021, 3:43pm UTC](https://discourse.julialang.org/t/target-variable-is-a-deterministic-function-of-random-variables/71935/6 "2021-11-25T15:43:00Z")
Yes, that’s correct. z = f(\text{random\_variable1, random\_variable2, ...}), with no error. If the values of the random variables are known, then the exact value of z is known.
Thanks for the suggestion, I’ve tried to use `Dirac` before and yes, the lack of proper sampling is an issue. Regarding `z ~ Normal(f(SomeRandomStuff), epsilon)`, are you suggesting I should remove the error term \epsilon from the equation y\_t = k \* x\_t + \epsilon, and instead add it into the z\_t equation z\_t = \text{Normal}(\text{lag\_transform}(\lambda, y\_t), \,\epsilon\_\text{new})? This could potentially work, however knowing the actual distribution of the original \epsilon is desirable.
Otherwise, if you meant to add an extra \epsilon term to the model, I’m not sure how this would be helpful. Is the idea to give it a very small prior, so that now the sampling works as intended and the added \epsilon has a negligible impact on the estimation of the other parameters?
---
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [November 25, 2021, 4:09pm UTC](https://discourse.julialang.org/t/target-variable-is-a-deterministic-function-of-random-variables/71935/7 "2021-11-25T16:09:27Z")
> [@Stephen2](#):
>
> Is the idea to give it a very small prior, so that now the sampling works as intended and the added \epsilon ϵ\epsilon has a negligible impact on the estimation of the other parameters?
Basically this, let’s say Z is typically O(1) then give epsilon a prior with scale say O(0.01) and let it loose. The sampled distribution will be not perfectly the one you want but it will sample. You can also take a large sample and then adjust the sample after. There are techniques for such things.
---
### Author:  [@Stephen2](https://discourse.julialang.org/u/Stephen2)
#### Post date: [November 29, 2021, 2:47pm UTC](https://discourse.julialang.org/t/target-variable-is-a-deterministic-function-of-random-variables/71935/8 "2021-11-29T14:47:43Z")
Excellent, this is producing very decent results! Are there any resources that go deeper into this? I would like to know more about how to choose an appropriate epsilon prior, as well as the other technique you mentioned.
For reference, this is my working model now:
```julia
@model function model_func(x, z)
k ~ Uniform(0.5, 3)
λ ~ Uniform(0.4, 1)
σ ~ truncated(Normal(0, 1), 0, Inf)
ϵ ~ Uniform(0.0001, 0.001)
y = [k * x_val + rand(Normal(0, σ)) for x_val in x]
powers = [λ ^ i for i in 1:5]
for n in 5:length(x)
z_val::Float32 =
powers[1] * y[n] + powers[2] * y[n - 1] + powers[3] * y[n - 2] + powers[4] * y[n - 3] + powers[5] * y[n - 4]
z[n-4] ~ Normal(z_val, ϵ)
end
end
```
---
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [November 29, 2021, 3:18pm UTC](https://discourse.julialang.org/t/target-variable-is-a-deterministic-function-of-random-variables/71935/9 "2021-11-29T15:18:01Z")
You could take a look at the CRC Handbook of Approximate Bayesian Computation (Sisson, Fan, and Beaumont).
Glad this is working for you!
From the perspective of Bayes as continuous valued logic, the Normal(z\_val,eps) likelihood is essentially saying that samples of z\_val which are very close to z are very credible values, and as you get farther away the credibility drops a lot. So you can think of it as an approximation to the model where they are exactly equal, or you can think of it as exactly a model in which z-z\_val differs by at most a very small amount. If you’re willing to take this last view, then there’s no adjustment required.
---
### Author:  [@Stephen2](https://discourse.julialang.org/u/Stephen2)
#### Post date: [November 29, 2021, 4:48pm UTC](https://discourse.julialang.org/t/target-variable-is-a-deterministic-function-of-random-variables/71935/10 "2021-11-29T16:48:29Z")
Great, thanks again!