# Sampling from the posterior with a constraint on the posterior mean

## Background

Under certain assumptions we know that being given the posterior mean and a family of conditional distributions, we can uniquely determine the joint distribution. I quote one of the theorems specifying appropriate conditions below. The theorem comes from this book and provides the conditions for power series distributions $$X|Y \sim PSD(Y)$$.

Theorem 7.2
Let $$(X,Y)$$ be a random vector such that either:
(a) $$S(X) = \{0,1,2,...,n\}$$ for some integer $$n$$ and the cardinality of $$S(Y)$$ is $$\leq n+2$$; or
(b) $$S(X)= \{0,1,2,...\}$$ and $$S(Y) \subseteq \{0,1,2,...\}$$.

Assume that the supports $$S(X)$$ and $$S(Y)$$ are known and that for any $$x \in S(X), y\in S(Y)$$ we have $$P(X=x|Y=y) = c(x)y^x/c^*(y),$$ where $$c$$ and $$c^*$$ are known. In addition, if $$S(Y)$$ is not bounded assume that $$\sum_{x \in S(X)} \sqrt[2x]{c(x)} = \infty.$$ Then the distribution of $$(X,Y)$$ is uniquely determined by $$\mathbb{E}(Y|X=x) = \psi(x), x\in S(x)$$.

## Example

Let $$X|Y=y \sim Pois(\lambda y)$$. Then $$X|Y \sim PSD(Y)$$ with $$c(x) = \lambda^x/x!$$ and $$c^*(y) = \exp(\lambda y)$$. $$\sum_{x=0}^{\infty} \sqrt[2x]{c(x)} \geq \sqrt{\lambda} \sum_{x=0}^{\infty} \frac{1}{x} = \infty$$, so the conditional expectation of $$Y$$ given $$X$$ uniquely determines the joint distribution.

We don't know anything about the joint distribution other than it exists and is unique. The proof is not constructive. I would like to find $$Y|X$$ by simulation methods.

In general let's assume we have the following setting - tractable likelihood $$X|\theta$$ and the functional form of the posterior mean $$\mathbb{E}(\theta|X) = \phi(X)$$.

## Question

Is it possible to simulate observations from the posterior and utilize the information about the posterior mean in an MCMC setting?

• "I would like to find 𝑌|𝑋 by simulation methods" How can you perform these simulations if you do not know the distribution of Y in your example? What are you gonna simulate? Or do you want to estimate Y|X given some sample of data? Apr 21 at 10:35
• @SextusEmpiricus Well, that's the question. I do not know the distribution of $Y$ but I know that it is unique as per the theorem and knowing the distribution of $X|Y$ and the functional form of $\mathbb{E}(Y|X) = \phi(X)$ is enough to find it. Apr 21 at 11:05

Given that both $$X$$ and $$Y$$ are integer valued, e.g. with a finite number of values, a likely trick in deriving the marginals of $$X$$ and $$Y$$ stands in writing enough linear relations between $$\mathbb P(X=x)$$ and $$\mathbb P(Y=y)$$ to identify the marginals. (After writing the derivation below, I checked the referenced book and found out that the authors reach the same conclusion in eqn (7.38).)
For instance, assume that $$S(X) =S(Y) = \{0,1,2,...,n\}$$. Then, for $$0\le \imath,\jmath\le n$$, $$\mathbb P(Y=\jmath|X=\imath) = \frac{\mathbb P(X=\imath|Y=\jmath)\mathbb P(Y=\jmath)}{\mathbb P(X=\imath)}=\frac{\jmath^\imath c(\imath)/c^*(\jmath)\overbrace{\mathbb P(Y=\jmath)}^{q(\jmath)}}{\underbrace{\mathbb P(X=\imath)}_{p(\imath)}}$$ Therefore, for $$0\le \imath\le n$$, $$\psi(\imath) = \sum_{\jmath=0}^n \jmath\mathbb P(Y=\jmath|X=\imath)=\sum_{\jmath=0}^n \jmath^{\imath+1} \frac{c(\imath)q(\jmath)}{c^*(\jmath)p(\imath)} = \frac{c(\imath)}{p(\imath)} \sum_{\jmath=0}^n \jmath^{\imath+1} \frac{q(\jmath)}{c^*(\jmath)}$$ or $$p(\imath)\frac{\psi(\imath)}{c(\imath)}=\sum_{\jmath=0}^n \jmath^{\imath+1} \frac{q(\jmath)}{c^*(\jmath)}\tag{1}$$ Furthermore, $$p(\imath)=\sum_{\jmath=0}^n \mathbb P(X=\imath|Y=\jmath) \mathbb P(Y=\jmath)=\sum_{\jmath=0}^n \frac{\jmath^\imath c(\imath)}{c^*(\jmath)}q(\jmath)\tag{2}$$ Hence, merging (1) and (2), $$\sum_{\jmath=0}^n \frac{\jmath^\imath c(\imath)}{c^*(\jmath)}q(j)\frac{\psi(\imath)}{c(\imath)}=\sum_{\jmath=0}^n \jmath^{i+1} \frac{q(\jmath)}{c^*(\jmath)}$$ or $$\sum_{\jmath=0}^n \frac{\jmath^\imath q(\jmath)\psi(\imath)}{c^*(\jmath)}=\sum_{\jmath=0}^n \jmath^{i+1} \frac{q(\jmath)}{c^*(\jmath)}$$ or again $$\sum_{\jmath=0}^n \left\{{\psi(\imath)}-{\jmath}\right\}\frac{\jmath^\imath}{c^*(\jmath)}q(\jmath)=0\tag{3}$$ which leads to a system of $$n+1$$ equations in $$\mathbf q=(q(0),\ldots,q(n))$$ and, along with the normalisation constraint $$\sum_{\jmath=0}^n q(\jmath)=1$$ it should lead to a unique derivation of the marginal distribution of $$Y$$.
Now, this does not directly answer the question about an MCMC resolution since solving (3) produces a marginal distribution and hence a way to directly simulate from the joint. (An unsubstantiated suggestion is to move $$\mathbf q$$ at each MCMC iteration by one gradient descent step when starting from a arbitrary value.)