# Turing.jl: fix some parameters to prior distribution

**URL:** <https://discourse.julialang.org/t/turing-jl-fix-some-parameters-to-prior-distribution/41095>\
**Category:** Probabilistic Programming\
**Created:** [June 9, 2020, 9:11pm UTC](https://discourse.julialang.org/t/turing-jl-fix-some-parameters-to-prior-distribution/41095 "2020-06-09T21:11:44Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![MFairley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mfairley/32/8599_2.png) [@MFairley](https://discourse.julialang.org/u/MFairley)\
**Post date:** [June 9, 2020, 9:11pm UTC](https://discourse.julialang.org/t/turing-jl-fix-some-parameters-to-prior-distribution/41095/1 "2020-06-09T21:11:44Z")

</div>

I’m testing out Turing.jl and I’m wondering if I can specify a subset of parameters to remain fixed at their prior distribution. I would still like to sample these parameters to get a joint distribution with the other parameters that I do want the posterior for.

For example, say I have the following model. Can I specify `β` to remain at its prior distribution?

```julia
@model bayes_sir(y) = begin
    # Calculate number of timepoints
    l = length(y)
    i₀ ~ Uniform(0.0,1.0)
    β ~ Uniform(0.0,1.0)
    I = i₀*1000.0
    u0=[1000.0-I,I,0.0,0.0]
    p=[β,10.0,0.25]
    tspan = (0.0,float(l))
    prob = ODEProblem(sir_ode!,
            u0,
            tspan,
            p)
    sol = solve(prob,
                Tsit5(),
                saveat = 1.0)
    sol_C = Array(sol)[4,:] # Cumulative cases
    sol_X = sol_C[2:end] - sol_C[1:(end-1)]
    l = length(y)
    for i in 1:l
      y[i] ~ Poisson(sol_X[i])
    end
end;

```

---

<div class="post-metadata">

**Author:** ![opera\_malenky](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/opera_malenky/32/8213_2.png) [@opera\_malenky](https://discourse.julialang.org/u/opera_malenky)\
**Post date:** [June 9, 2020, 9:24pm UTC](https://discourse.julialang.org/t/turing-jl-fix-some-parameters-to-prior-distribution/41095/2 "2020-06-09T21:24:08Z")

</div>

I don’t know the answer, but you mean something like the [cut function](http://statistik.wu-wien.ac.at/research/talks/resources/Plummer_WU_2015.pdf) in WinBUGS?

---

<div class="post-metadata">

**Author:** ![MFairley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mfairley/32/8599_2.png) [@MFairley](https://discourse.julialang.org/u/MFairley)\
**Post date:** [June 9, 2020, 9:33pm UTC](https://discourse.julialang.org/t/turing-jl-fix-some-parameters-to-prior-distribution/41095/3 "2020-06-09T21:33:08Z")

</div>

That looks similar but perhaps my goal is simpler. I don’t want to update the distribution of a subsets of parameters at all, not even using a subset of the data.

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [June 9, 2020, 9:45pm UTC](https://discourse.julialang.org/t/turing-jl-fix-some-parameters-to-prior-distribution/41095/4 "2020-06-09T21:45:57Z")

</div>

Every ~ indicates that we perform inference for that parameter and that we account for it in the joint probability. But you can simply draw beta using the rand function ~~and add its log probability to the joint.~~

Does that help?  
May I also ask why you don’t want to update the prior?

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [June 9, 2020, 9:56pm UTC](https://discourse.julialang.org/t/turing-jl-fix-some-parameters-to-prior-distribution/41095/5 "2020-06-09T21:56:41Z")

</div>

You can use the approach outlined here [Is there a turing alternative to pm.Deterministic from pymc3? - #3 by torfjelde](https://discourse.julialang.org/t/is-there-a-turing-alternative-to-pm-deterministic-from-pymc3/38667/3) to keep track of beta.

---

<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:** [June 9, 2020, 9:58pm UTC](https://discourse.julialang.org/t/turing-jl-fix-some-parameters-to-prior-distribution/41095/6 "2020-06-09T21:58:05Z")

</div>

You’re trying to do inference on a **subset** of the parameters? Like if you have parameters a,b you want to keep the distribution of a as the prior, but do inference on b?

You could do this by sampling from the “fixed” parameters, then for each of those fixed sample values, run inference on b conditional on a remaining equal to the sample value… When you pool all the runs, the marginal distribution of a would still be the prior, but the b would be according to p(b | a, data) for each sample value of a.

---

<div class="post-metadata">

**Author:** ![MFairley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mfairley/32/8599_2.png) [@MFairley](https://discourse.julialang.org/u/MFairley)\
**Post date:** [June 9, 2020, 10:07pm UTC](https://discourse.julialang.org/t/turing-jl-fix-some-parameters-to-prior-distribution/41095/7 "2020-06-09T22:07:52Z")

</div>

> [@trappmartin](#):
>
> Every ~ indicates that we perform inference for that parameter and that we account for it in the joint probability. But you can simply draw beta using the rand function ~~and add its log probability to the joint.~~
> 
> Does that help?  
> May I also ask why you don’t want to update the prior?

Great, so something like this:

```julia
@model bayes_sir(y) = begin
    # Calculate number of timepoints
    l = length(y)
    i₀ ~ Uniform(0.0,1.0)
    β = rand(Uniform(0.0,1.0))
    I = i₀*1000.0
    u0=[1000.0-I,I,0.0,0.0]
    p=[β,10.0,0.25]
    tspan = (0.0,float(l))
    prob = ODEProblem(sir_ode!,
            u0,
            tspan,
            p)
    sol = solve(prob,
                Tsit5(),
                saveat = 1.0)
    sol_C = Array(sol)[4,:] # Cumulative cases
    sol_X = sol_C[2:end] - sol_C[1:(end-1)]
    l = length(y)
    for i in 1:l
      y[i] ~ Poisson(sol_X[i])
    end
end;

```

and then use the trick from [Is there a turing alternative to pm.Deterministic from pymc3? - #3 by torfjelde](https://discourse.julialang.org/t/is-there-a-turing-alternative-to-pm-deterministic-from-pymc3/38667/3) to keep track of the samples from beta?

There are a couple of applications I had in mind:

1. For computational efficiency, when I know that a parameter has no influence on the likelihood.

2. When the prior was estimated using some external source that I think is more reliable than my own data. For example, say I observe deaths in patients from treatments A and B. A clinical trial reports a hazard rate ratio (hr) for the effectiveness of A vs. B. I assume the likelihood is: deaths\_a ~ Poisson(lambda \* hr) and deaths\_b ~ Poisson(lambda). I want to run inference on lambda but I want to leave hr unchanged from its prior.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [June 9, 2020, 10:23pm UTC](https://discourse.julialang.org/t/turing-jl-fix-some-parameters-to-prior-distribution/41095/8 "2020-06-09T22:23:13Z")

</div>

I think we should be able to formally support this use case by allowing the `Prior` inference algorithm to be a `GibbsComponent`. Then you can use `Gibbs` to specify the splitting.

---

<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:** [June 9, 2020, 10:43pm UTC](https://discourse.julialang.org/t/turing-jl-fix-some-parameters-to-prior-distribution/41095/9 "2020-06-09T22:43:13Z")

</div>

> [@MFairley](#):
>
> I want to run inference on lambda but I want to leave hr unchanged from its prior.

I don’t think that’s a good example. What you’d do here is use the prior information from the previous reports, to provide a very tight prior, and then you’d include the `hr` in your parameter set and do inference on it. If there’s not much information in your data, it won’t change much from the prior, on the other hand if there is information in your data set, it will change from the prior, and that’s appropriate. Basically a tight prior from a previous study is still just a prior not a set-in-stone fact (unless the value is _very_ tight, in which case you might as well just enter a number directly into your model)

---

<div class="post-metadata">

**Author:** ![trappmartin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trappmartin/32/1165_2.png) [@trappmartin](https://discourse.julialang.org/u/trappmartin)\
**Post date:** [June 10, 2020, 7:44am UTC](https://discourse.julialang.org/t/turing-jl-fix-some-parameters-to-prior-distribution/41095/10 "2020-06-10T07:44:49Z")

</div>

I agree with @dlakelan. If you have strong prior assumptions, it is better to incorporate these into the generative process using a tight prior or logical constraints.

---

<div class="post-metadata">

**Author:** ![MFairley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mfairley/32/8599_2.png) [@MFairley](https://discourse.julialang.org/u/MFairley)\
**Post date:** [June 10, 2020, 3:35pm UTC](https://discourse.julialang.org/t/turing-jl-fix-some-parameters-to-prior-distribution/41095/11 "2020-06-10T15:35:57Z")

</div>

> [@trappmartin](#):
>
> I agree with @dlakelan. If you have strong prior assumptions, it is better to incorporate these into the generative process using a tight prior or logical constraints.

There could be a wide prior for `hr` estimated from the clinical trial but the clinical trial is an unbiased estimate of `hr` because patients were randomized to treatment arms. My own, observational data may have some unknown selection effects so I want to build a model that constrains `hr` to be consistent with the results of the clinical trial, and not let my own data influence `hr`. I think this idea is similar to the reasons discussed here [http://statistik.wu-wien.ac.at/research/talks/resources/Plummer\_WU\_2015.pdf](http://statistik.wu-wien.ac.at/research/talks/resources/Plummer_WU_2015.pdf) except that the data `Z` is not available and I only have access to the distribution for `φ` already estimated using `Z` by somebody else.

---

<div class="post-metadata">

**Author:** ![MFairley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mfairley/32/8599_2.png) [@MFairley](https://discourse.julialang.org/u/MFairley)\
**Post date:** [June 10, 2020, 4:32pm UTC](https://discourse.julialang.org/t/turing-jl-fix-some-parameters-to-prior-distribution/41095/12 "2020-06-10T16:32:33Z")

</div>

> [@trappmartin](#):
>
> Every ~ indicates that we perform inference for that parameter and that we account for it in the joint probability. But you can simply draw beta using the rand function

I don’t think this approach will work because inference on beta still occurs. For example:

```julia
@model bayes_sir(y) = begin
    # Calculate number of timepoints
    l = length(y)
    i₀ ~ Uniform(0.0,1.0)
    # β ~ Uniform(0.0,1.0)
    β ~ Determin(rand(Uniform(0.0, 1.0)))
    I = i₀*1000.0
    u0=[1000.0-I,I,0.0,0.0]
    p=[β,10.0,0.25]
    tspan = (0.0,float(l))
    prob = ODEProblem(sir_ode!,
            u0,
            tspan,
            p)
    sol = solve(prob,
                Tsit5(),
                saveat = 1.0)
    sol_C = Array(sol)[4,:] # Cumulative cases
    sol_X = sol_C[2:end] - sol_C[1:(end-1)]
    l = length(y)
    for i in 1:l
      y[i] ~ Poisson(sol_X[i])
    end
end;

ode_smc = sample(bayes_sir(generate_data(1234)), SMC(), 10000);

```

returns

```julia
Quantiles
  parameters 2.5% 25.0% 50.0% 75.0% 97.5%
  ────────── ────── ────── ────── ────── ──────
          i₀ 0.0107 0.0107 0.0107 0.0192 0.0192
           β 0.0435 0.0435 0.0501 0.0501 0.0514

```

So beta is no longer at its prior distribution. I guess I would need to do something like the 2-stage process @dlakelan described:

> [@dlakelan](#):
>
> You could do this by sampling from the “fixed” parameters, then for each of those fixed sample values, run inference on b conditional on a remaining equal to the sample value… When you pool all the runs, the marginal distribution of a would still be the prior, but the b would be according to p(b | a, data) for each sample value of a.

I imagine this would be very computationally expensive if I implemented myself because I would have to re-run the Turing sampler for each sample of beta.

---

<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:** [June 10, 2020, 4:38pm UTC](https://discourse.julialang.org/t/turing-jl-fix-some-parameters-to-prior-distribution/41095/13 "2020-06-10T16:38:55Z")

</div>

> [@MFairley](#):
>
> There could be a wide prior for `hr` estimated from the clinical trial but the clinical trial is an unbiased estimate of `hr` because patients were randomized to treatment arms. My own, observational data may have some unknown selection effects so I want to build a model that constrains `hr` to be consistent with the results of the clinical trial, and not let my own data influence `hr`

I think the issue is that you’re building the wrong kind of model.

If you think the clinical trial gave a good unbiased estimate of a distribution for beta, and that your data has a bias… Then include a parameter in your model that represents the bias. Instead of using beta to describe your data, use beta + bias.

Another way of saying this is “partial pooling” or “hierarchical modeling”.

beta ~ my\_clinical\_prior  
mystudybeta ~ my\_bias\_size(beta)

mydata ~ my\_likelihood(mystudybeta)

This is entirely the same concept, except it is also easier to incorporate logical constraints if you formulate it this way.

---

<div class="post-metadata">

**Author:** ![MFairley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mfairley/32/8599_2.png) [@MFairley](https://discourse.julialang.org/u/MFairley)\
**Post date:** [June 10, 2020, 5:20pm UTC](https://discourse.julialang.org/t/turing-jl-fix-some-parameters-to-prior-distribution/41095/14 "2020-06-10T17:20:06Z")

</div>

Okay, but I’m still interested to know how I could write code to force beta to stay at its prior.

I guess as @mohamed82008 said, this could be achieved by something like:

```julia
chn = sample(bayes_sir(generate_data(1234)), Gibbs(NUTS(-1, 0.65, :i₀), Prior(:β)), 1000)

```

---

<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:** [June 10, 2020, 5:27pm UTC](https://discourse.julialang.org/t/turing-jl-fix-some-parameters-to-prior-distribution/41095/15 "2020-06-10T17:27:35Z")

</div>

Think about what it means for Beta to stay at its prior distribution with respect to proper Bayesian inference. It means that the likelihood function is a constant with respect to Beta. And this in turn means that Beta doesn’t inform you about what your data might be.

All you have to do is simply don’t use Beta 😉

As the distribution for the biased beta tends towards infinitely wide… The distribution for the unbiased Beta tends towards the prior.

There is no way to do actual Bayesian inference where you use the distribution for Beta and you don’t distort it in the posterior. If you make the bias infinitely wide, you don’t use the Beta, if you make the bias less than infinitely wide you use the Beta information, but you distort the Beta somewhat. It’s a continuum of possibilities. Simply choose that place on the continuum where you’re happy.

Think about it like this… you have three balls connected by springs. The ball on the right is your data, C, the ball in the middle is your biased parameter, B, and the ball on the left is your clinical Beta, A. Like this:

A ---- B ----- C

If you make the spring between A,B very very soft… then if you pull on C to put it where the real data lies… you’ll bring B along with you, but A will stay behind because B isn’t pulling on it very much…

If you tighten up the spring between A, B you’ll pull A along with B somewhat. If you make A,B connected by a rigid rod with zero length, A will come along and always be exactly where B is.

The model you started with is where A,B are connected by a rigid rod with zero length… By partially decoupling A,B using a “soft spring” (namely a prior on B that is wide and centered around the location of A) then A will remain mainly where it was. If you completely decouple by providing some other prior for B that doesn’t involve A at all… then you’ll sample from the clinical prior for A… but it won’t affect your choice of B at all.