# Computing or sampling a posterior with samples observed through a dimensional reduction transformation

**URL:** https://discourse.mc-stan.org/t/computing-or-sampling-a-posterior-with-samples-observed-through-a-dimensional-reduction-transformation/36257
**Category:** Modeling
**Created:** [August 15, 2024, 4:47am UTC](https://discourse.mc-stan.org/t/computing-or-sampling-a-posterior-with-samples-observed-through-a-dimensional-reduction-transformation/36257 "2024-08-15T04:47:12Z")
**Posts on this page:** 2
**Page:** 1

<div class="post-metadata">

### Author: ![Latiao](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/latiao/32/18989_2.png) [@Latiao](https://discourse.mc-stan.org/u/Latiao)
#### Post date: [August 15, 2024, 4:47am UTC](https://discourse.mc-stan.org/t/computing-or-sampling-a-posterior-with-samples-observed-through-a-dimensional-reduction-transformation/36257/1 "2024-08-15T04:47:12Z")

</div>

Let \boldsymbol \theta be a vector of parameters, with a known prior \pi(\boldsymbol \theta).  
Let \boldsymbol x\_1,...,\boldsymbol x\_n be i.i.d. samples with \boldsymbol x|\boldsymbol \theta known.  
The problem is that we do not observe \boldsymbol x\_i, but \boldsymbol g(\boldsymbol x\_i). Here \boldsymbol g is a dimensional reduction transformation. How to numerically compute or sample the posterior \pi(\boldsymbol \theta | \boldsymbol g(\boldsymbol x\_1),...,g(\boldsymbol x\_n)) efficiently?

A simple example is here. For a rectangle, the width w\sim {\rm Uniform}(0, a) and the length l\sim {\rm Uniform}(0,b), and (a,b) follows a known prior.  
But we cannot observe (w,l). Instead, we can observe the area s = w\cdot l.  
How to numerically compute or sample the posterior \pi(a,b|s)?

Related problem: [Estimating a posterior for a parameter which is only observed through (a part of a) linear transformation](https://discourse.mc-stan.org/t/estimating-a-posterior-for-a-parameter-which-is-only-observed-through-a-part-of-a-linear-transformation/6309)

---

<div class="post-metadata">

### Author: ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)
#### Post date: [August 22, 2024, 7:04am UTC](https://discourse.mc-stan.org/t/computing-or-sampling-a-posterior-with-samples-observed-through-a-dimensional-reduction-transformation/36257/2 "2024-08-22T07:04:14Z")

</div>

Do you assume that you then observe \boldsymbol{g}(\boldsymbol{x}\_i) without any additional noise?

If there is noise in observing \boldsymbol{g}(\boldsymbol{x}\_i) that - at least in principle - you can treat \boldsymbol x\_1,...,\boldsymbol x\_n as explicit parameters in your model and directly include the dimensionality reduction in your model to compute a the distribution \pi(\boldsymbol{g}(\boldsymbol{x}\_i) | \boldsymbol{x}\_i)

If there is no noise, then you need to somehow invert the dimensional reduction transformation, i.e. get a description of the set \{\boldsymbol x\_1,...,\boldsymbol x\_n | \boldsymbol{g} (\boldsymbol x\_1) = y\_1,...,\boldsymbol{g}(\boldsymbol x\_n) = y\_n\} - this might be possible for simple methods like using first few PCA components, but will be pretty annoying.

E.g. in the area example you would treat w as an extra unknown parameter and compute l = \frac{s}{w}. You can now directly express your model for \pi(w,l | \theta) while adding a Jacobian adjustment for the transform.

An approximate approach might be possible by doing something similar to multiple imputation (since your problem can be understood as a missing data problem). I.e. you don’t insist on finding all \boldsymbol x that produce the observed reduced data, but just get a bunch of different possible values, and for each of those fit a model for \pi(\boldsymbol x | \boldsymbol \theta) and then pool the samples from all of those fits.

You could likely obtain suitable samples of \boldsymbol x\_1,...,\boldsymbol x\_n by starting with random values for \boldsymbol {x}^\* and then run an optimization algorithm to minimize \sum\_i^n (\boldsymbol{g}(\boldsymbol x\_n) - \boldsymbol{g}(\boldsymbol {x}^\*\_n))^2 or similar.

In all cases I think the problem is likely to be computationally somewhat ill-posed - there might be many different regions of \boldsymbol \theta that produce similarly looking distribution over \boldsymbol g(x) and thus you will have multimodality.

Hope that clarifies more than confuses and best of luck with your model!
