# Implementing model where data is a function of random variables

**URL:** <https://discourse.julialang.org/t/implementing-model-where-data-is-a-function-of-random-variables/75354>\
**Category:** Statistics\
**Tags:** turing\
**Created:** [January 28, 2022, 11:34am UTC](https://discourse.julialang.org/t/implementing-model-where-data-is-a-function-of-random-variables/75354 "2022-01-28T11:34:50Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Pablo\_Marchant](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pablo_marchant/32/32828_2.png) [@Pablo\_Marchant](https://discourse.julialang.org/u/Pablo_Marchant)\
**Post date:** [January 28, 2022, 11:34am UTC](https://discourse.julialang.org/t/implementing-model-where-data-is-a-function-of-random-variables/75354/1 "2022-01-28T11:34:50Z")

</div>

Hi,

I’ve got a question regarding implementation of a model where an observed variable can be described as a function of other random variables. My issues here are perhaps just due to my inexperience in using probabilistic programming rather than a problem specific to Turing, so any help is appreciated. In particular I have a problem where an observed quantity I is given by a function of the form:

I=A\cos(\phi)+B

where \phi\sim \mathcal{N}(\mu\_1,\sigma\_1) and \phi\sim \mathcal{N}(\mu\_2,\sigma\_2). Model parameters are then A,\mu\_1,\sigma\_1,\mu\_2,\sigma\_2. I could in principle derive the distribution for I, I\sim\mathcal{F}(A,\mu\_1,\sigma\_1,\mu\_2,\sigma\_2), and construct a Turing model like

```julia
@model function offset_model(y)
	A ~ Uniform(0,100)
	mu1 ~ Uniform(0,100)
	sigma1 ~ Uniform(0,100)
	mu2 ~ Uniform(0,100)
	sigma2 ~ Uniform(0,100)
    # The number of observations.
    N = length(y)
    for n in 1:N
        y[n] ~ F(A,mu1,sigma1,mu2,sigma2)
    end
end;

```

And using this, for a set of observations of I, I could then draw the posteriors for my model parameters.

But I’m wondering if there is some way to do this automatically in julia+turing without having to derive what would be a pretty complex distribution function for I. In particular, a flexible way to specify things would be useful as \phi or B can potentially follow different distributions, and having to derive the corresponding distributions for I can be very cumbersome.

As said before, any help is greatly appreciated!

---

<div class="post-metadata">

**Author:** ![blackeneth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/blackeneth/32/10353_2.png) [@blackeneth](https://discourse.julialang.org/u/blackeneth)\
**Post date:** [January 30, 2022, 8:43pm UTC](https://discourse.julialang.org/t/implementing-model-where-data-is-a-function-of-random-variables/75354/2 "2022-01-30T20:43:03Z")

</div>

You could possibly use some functionality from [MixedModels.jl Documentation · MixedModels](https://juliastats.org/MixedModels.jl/stable/)

---

<div class="post-metadata">

**Author:** ![sethaxen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sethaxen/32/35604_2.png) [@sethaxen](https://discourse.julialang.org/u/sethaxen)\
**Post date:** [January 31, 2022, 4:52pm UTC](https://discourse.julialang.org/t/implementing-model-where-data-is-a-function-of-random-variables/75354/3 "2022-01-31T16:52:43Z")

</div>

> [@Pablo\_Marchant](#):
>
> In particular I have a problem where an observed quantity I II is given by a function of the form:
> 
> I=A\cos(\phi)+B I=Acos(ϕ)+BI=A\cos(\phi)+B
> 
> where \phi\sim \mathcal{N}(\mu\_1,\sigma\_1) ϕ∼N(μ1,σ1)\phi\sim \mathcal{N}(\mu\_1,\sigma\_1) and \phi\sim \mathcal{N}(\mu\_2,\sigma\_2) ϕ∼N(μ2,σ2)\phi\sim \mathcal{N}(\mu\_2,\sigma\_2) .

Is there a notational error? Or, if not, what does it mean that \phi is drawn from two different distributions here?

---

<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:** [January 31, 2022, 5:23pm UTC](https://discourse.julialang.org/t/implementing-model-where-data-is-a-function-of-random-variables/75354/4 "2022-01-31T17:23:50Z")

</div>

> [@Pablo\_Marchant](#):
>
> In particular I have a problem where an observed quantity I II is given by a function of the form:
> 
> I=A\cos(\phi)+B I=Acos(ϕ)+B

if you allow a small error term, then you can try to do something like this:

```julia
y[i] ~ normal(A*cos(phi)+B, small_error)

```

I suppose maybe I is typicall O(1) or some such thing, so then let small\_error be something like .001 to get started and see what happens.

---

<div class="post-metadata">

**Author:** ![Pablo\_Marchant](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pablo_marchant/32/32828_2.png) [@Pablo\_Marchant](https://discourse.julialang.org/u/Pablo_Marchant)\
**Post date:** [January 31, 2022, 5:27pm UTC](https://discourse.julialang.org/t/implementing-model-where-data-is-a-function-of-random-variables/75354/5 "2022-01-31T17:27:31Z")

</div>

Oh, yes, that was dumb of me. I have that \phi\sim \mathcal{N}(\mu\_1,\sigma\_1) and B\sim \mathcal{N}(\mu\_2,\sigma\_2) (so the second normal is for B, not \phi).

---

<div class="post-metadata">

**Author:** ![Pablo\_Marchant](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pablo_marchant/32/32828_2.png) [@Pablo\_Marchant](https://discourse.julialang.org/u/Pablo_Marchant)\
**Post date:** [January 31, 2022, 8:00pm UTC](https://discourse.julialang.org/t/implementing-model-where-data-is-a-function-of-random-variables/75354/6 "2022-01-31T20:00:22Z")

</div>

Did not know about MixedModels, thanks for pointing it out! This actually helps me realize that what I’m trying to do is to construct a non-linear mixed effects model, which upon a bit of inspection turns out to be a much more complex problem than I thought at first.