Simulating $f(Ax + b)$ Let $f_1$ be the density function of the logistic distribution, and let the product density function $f(x,y,z) = f_1(x)f_1(y)f_1(z)$. A $3 \times 2$-matrix $A$ is given. Its two columns are orthogonal and each of them is normalized. A vector $b$ of length $3$ is given as well.
I want to simulate the distribution on $\mathbb{R}^2$ whose density function is, up to a normalization constant, the function $(u,v) \mapsto f\bigl(A\cdot(u,v) + b\bigr)$. What are the possible methods?
 A: Let $Y = A^{\dagger}(X-b)$ where $A^{\dagger}$ is the pseudo-inverse of $A$. If $X\sim f(x)$, then we have $Y\sim f(Ay+b)$
Proof: let $g:\mathbb{R}^d \to \mathbb{R}$ be a function with finite support,
$$\mathbb{E}[g(Y)]=\int_{\mathbb{R}^d} g(A^{\dagger}(x-b)) f(x)dx=\int_{\mathbb{R}^d} g(A^{-1}(Ay+b-b)) f(Ay+b)d(Ay+b)$$
Hence (this is a tricky part),
$$\mathbb{E}[g(Y)]=\int_{\mathbb{R}^d} g(y) f(Ay+b)|\mathrm{det}(A)|dy$$
where $\mathrm{det}(A)=\prod \alpha_i$ where $\alpha_i$ are the non-zero singular values of $A$ (using for example 3.6.2 in Measure Theory, Vol 1, by V.I. Bogachev, I think this work I am not sure, I did not do the whole computation).
Illustration:

code:
import numpy as np
from scipy.linalg import pinv

Sigma = np.diag([1, 2, 3])
X = np.random.multivariate_normal([0,0,0], Sigma, size=1000)

A = np.array([[1,2], [3,4], [1,1]])
b = np.array([1,2,3])

Y = np.array([pinv(A)@(x-b) for x in X])

import matplotlib.pyplot as plt
f, (ax1, ax2) = plt.subplots(1, 2, sharex=True, sharey=True)

ax1.hist2d(Y[:,0], Y[:,1], bins=(50, 50), cmap=plt.cm.Reds)

from scipy.stats import multivariate_normal
law = multivariate_normal([0,0,0], Sigma)
xs = np.linspace(-10, 6, 50)
ys = np.linspace(-5, 7, 50)
z = np.array([ law.pdf(A@np.array([x,y]) +b ) for y in ys for x in xs])

X, Y = np.meshgrid(xs, ys)
Z = z.reshape(50, 50)
ax2.pcolor(X, Y, Z, cmap=plt.cm.Reds)
plt.show()

