|
| 1 | +--- |
| 2 | +title: Tracking Extra Quantities |
| 3 | +engine: julia |
| 4 | +aliases: |
| 5 | + - ../../tutorials/usage-generated-quantities/index.html |
| 6 | + - ../generated-quantities/index.html |
| 7 | +--- |
| 8 | + |
| 9 | +```{julia} |
| 10 | +#| echo: false |
| 11 | +#| output: false |
| 12 | +using Pkg; |
| 13 | +Pkg.instantiate(); |
| 14 | +``` |
| 15 | + |
| 16 | +Often, there are quantities in models that we might be interested in viewing the values of, but which are not random variables in the model that are explicitly drawn from a distribution. |
| 17 | + |
| 18 | +As a motivating example, the most natural parameterization for a model might not be the most computationally feasible. |
| 19 | +Consider the following (efficiently reparametrized) implementation of Neal's funnel [(Neal, 2003)](https://arxiv.org/abs/physics/0009028): |
| 20 | + |
| 21 | +```{julia} |
| 22 | +using Turing |
| 23 | +setprogress!(false) |
| 24 | +
|
| 25 | +@model function Neal() |
| 26 | + # Raw draws |
| 27 | + y_raw ~ Normal(0, 1) |
| 28 | + x_raw ~ arraydist([Normal(0, 1) for i in 1:9]) |
| 29 | +
|
| 30 | + # Transform: |
| 31 | + y = 3 * y_raw |
| 32 | + x = exp.(y ./ 2) .* x_raw |
| 33 | + return nothing |
| 34 | +end |
| 35 | +``` |
| 36 | + |
| 37 | +In this case, the random variables exposed in the chain (`x_raw`, `y_raw`) are not in a helpful form — what we're after are the deterministically transformed variables `x` and `y`. |
| 38 | + |
| 39 | +There are two ways to track these extra quantities in Turing.jl. |
| 40 | + |
| 41 | +## Using `:=` (during inference) |
| 42 | + |
| 43 | +The first way is to use the `:=` operator, which behaves exactly like `=` except that the values of the variables on its left-hand side are automatically added to the chain returned by the sampler. |
| 44 | +For example: |
| 45 | + |
| 46 | +```{julia} |
| 47 | +@model function Neal_coloneq() |
| 48 | + # Raw draws |
| 49 | + y_raw ~ Normal(0, 1) |
| 50 | + x_raw ~ arraydist([Normal(0, 1) for i in 1:9]) |
| 51 | +
|
| 52 | + # Transform: |
| 53 | + y := 3 * y_raw |
| 54 | + x := exp.(y ./ 2) .* x_raw |
| 55 | +end |
| 56 | +
|
| 57 | +sample(Neal_coloneq(), NUTS(), 1000) |
| 58 | +``` |
| 59 | + |
| 60 | +## Using `returned` (post-inference) |
| 61 | + |
| 62 | +Alternatively, one can specify the extra quantities as part of the model function's return statement: |
| 63 | + |
| 64 | +```{julia} |
| 65 | +@model function Neal_return() |
| 66 | + # Raw draws |
| 67 | + y_raw ~ Normal(0, 1) |
| 68 | + x_raw ~ arraydist([Normal(0, 1) for i in 1:9]) |
| 69 | +
|
| 70 | + # Transform and return as a NamedTuple |
| 71 | + y = 3 * y_raw |
| 72 | + x = exp.(y ./ 2) .* x_raw |
| 73 | + return (x=x, y=y) |
| 74 | +end |
| 75 | +
|
| 76 | +chain = sample(Neal_return(), NUTS(), 1000) |
| 77 | +``` |
| 78 | + |
| 79 | +The sampled chain does not contain `x` and `y`, but we can extract the values using the `returned` function. |
| 80 | +Calling this function outputs an array: |
| 81 | + |
| 82 | +```{julia} |
| 83 | +nts = returned(Neal_return(), chain) |
| 84 | +``` |
| 85 | + |
| 86 | +where each element of which is a NamedTuple, as specified in the return statement of the model. |
| 87 | + |
| 88 | +```{julia} |
| 89 | +nts[1] |
| 90 | +``` |
| 91 | + |
| 92 | +## Which to use? |
| 93 | + |
| 94 | +There are some pros and cons of using `returned`, as opposed to `:=`. |
| 95 | + |
| 96 | +Firstly, `returned` is more flexible, as it allows you to track any type of object; `:=` only works with variables that can be inserted into an `MCMCChains.Chains` object. |
| 97 | +(Notice that `x` is a vector, and in the first case where we used `:=`, reconstructing the vector value of `x` can also be rather annoying as the chain stores each individual element of `x` separately.) |
| 98 | + |
| 99 | +A drawback is that naively using `returned` can lead to unnecessary computation during inference. |
| 100 | +This is because during the sampling process, the return values are also calculated (since they are part of the model function), but then thrown away. |
| 101 | +So, if the extra quantities are expensive to compute, this can be a problem. |
| 102 | + |
| 103 | +To avoid this, you will essentially have to create two different models, one for inference and one for post-inference. |
| 104 | +The simplest way of doing this is to add a parameter to the model argument: |
| 105 | + |
| 106 | +```{julia} |
| 107 | +@model function Neal_coloneq_optional(track::Bool) |
| 108 | + # Raw draws |
| 109 | + y_raw ~ Normal(0, 1) |
| 110 | + x_raw ~ arraydist([Normal(0, 1) for i in 1:9]) |
| 111 | +
|
| 112 | + if track |
| 113 | + y = 3 * y_raw |
| 114 | + x = exp.(y ./ 2) .* x_raw |
| 115 | + return (x=x, y=y) |
| 116 | + else |
| 117 | + return nothing |
| 118 | + end |
| 119 | +end |
| 120 | +
|
| 121 | +chain = sample(Neal_coloneq_optional(false), NUTS(), 1000) |
| 122 | +``` |
| 123 | + |
| 124 | +The above ensures that `x` and `y` are not calculated during inference, but allows us to still use `returned` to extract them: |
| 125 | + |
| 126 | +```{julia} |
| 127 | +returned(Neal_coloneq_optional(true), chain) |
| 128 | +``` |
| 129 | + |
| 130 | +Another equivalent option is to use a submodel: |
| 131 | + |
| 132 | +```{julia} |
| 133 | +@model function Neal() |
| 134 | + y_raw ~ Normal(0, 1) |
| 135 | + x_raw ~ arraydist([Normal(0, 1) for i in 1:9]) |
| 136 | + return (x_raw=x_raw, y_raw=y_raw) |
| 137 | +end |
| 138 | +
|
| 139 | +chain = sample(Neal(), NUTS(), 1000) |
| 140 | +
|
| 141 | +@model function Neal_with_extras() |
| 142 | + neal ~ to_submodel(Neal(), false) |
| 143 | + y = 3 * neal.y_raw |
| 144 | + x = exp.(y ./ 2) .* neal.x_raw |
| 145 | + return (x=x, y=y) |
| 146 | +end |
| 147 | +
|
| 148 | +returned(Neal_with_extras(), chain) |
| 149 | +``` |
| 150 | + |
| 151 | +Note that for the `returned` call to work, the `Neal_with_extras()` model must have the same variable names as stored in `chain`. |
| 152 | +This means the submodel `Neal()` must not be prefixed, i.e. `to_submodel()` must be passed a second parameter `false`. |
0 commit comments