Remark:
I've added a detailed description of the actual setting of the application to the end of the question.
I'm running the Metropolis-Hastings algorithm with target distribution $\hat\mu$ (see definitions below) and proposal kernel $\hat Q$ on the product state space $\hat E:=I\times E$ and need to estimate the asymptotic variance of an estimator build up on the generated chain $(\hat X_n)_{n\in\mathbb N_0}=((S_n,X_n))_{n\in\mathbb N_0}$ and the corresponding proposal sequence $(\hat Y_n)_{n\in\mathbb N}=((T_n,Y_n))_{n\in\mathbb N}$.
To be precise, let \begin{equation}\begin{split}\hat\rho((i,x),(j,y))&:=\left.\begin{cases}\displaystyle\frac{\hat p(j,y)}{\hat q((i,x),(j,y))}&\text{, if }w_j(y)>0\\0&\text{, otherwise}\end{cases}\right\}\\&=\left.\begin{cases}\displaystyle\frac{p(y)}{q_i(y)u'\left(\varphi_i^{-1}(x),\varphi_i^{-1}(y)\right)}&\text{, if }w_j(y)>0\\0&\text{, otherwise}\end{cases}\right\}\end{split}\end{equation} for $(i,x),(j,y)\in\hat E$ and note that $$I_ng:=\frac{\sum_{i=1}^n\hat\rho\left(\hat X_{i-1},\hat Y_i\right)g(Y_i)}{\sum_{i=1}^n\hat\rho\left(\hat X_{i-1},\hat Y_i\right)}\;\;\;\text{for }n\in\mathbb N$$ is an asymptotically consistent estimator for $\mu g$, $g\in\mathcal L^1(\mu)$. If $g\in\mathcal L^1(\mu)$ and $g_0:=g-\mu g$ satisfy $$\sigma^2(g):=\int\hat\mu({\rm d}(i,x))\int\mu({\rm d}y)\hat\rho((i,x),y)|g_0(y)|^2<\infty\tag5,$$ then $$\sqrt nI_ng_0\xrightarrow{n\to\infty}\mathcal N(0,\sigma^2(g))\tag6.$$
Now, for a fixed $\mathcal E$-measurable $g:E\to[0,\infty)$ with $\{p=0\}\subseteq\{g=0\}$, $\lambda g<\infty$, $$(\iota g)(x):=\left.\begin{cases}\displaystyle\frac{g(x)}{p(x)}&\text{, if }p(x)>0\\0&\text{, otherwise}\end{cases}\right\}\;\;\;\text{for }x\in E$$ and $g_0:=g-\lambda g$, I need to do the following: First note that \begin{equation}\begin{split}\sigma^2(\iota g)&=\sum_{i\in I}\int\mu({\rm d}x)w_i(x)\underbrace{\int_{\{\:p\:>\:0\:\}}\lambda({\rm d}y)\frac{|g_0(y)|^2}{q_i(y)u'\left(\varphi_i^{-1}(x),\varphi_i^{-1}(y)\right)}}_{=:\:a_i(x)}\\&=\sum_{i\in I}\langle a_i,w_i\rangle_{L^2(\mu)}.\end{split}\tag7\end{equation} From the right-hand side of $(7)$ it is immediate that a choice of $(w_i)_{i\in I}$ minimizing $\sigma^2(\iota g)$ is given by $$w_i(x):=\left.\begin{cases}1&\text{, if }i=i_\ast(x)\\0&\text{, otherwise}\end{cases}\right\}\;\;\;\text{for }(i,x)\in E\hat,\tag0$$ where $$i_\ast(x):=\min\underset{i\in I}{\operatorname{arg min}}a_i(x)\;\;\;\text{for }x\in E$$ ($\operatorname{arg min}$ is considered as being set-valued).
Question: How can I compute $w_i(x)$, $x\in E$, numerically?
EDIT: Actually, I don't necessarily need to compute $w_i(x)$, $x\in E$, directly, but to sample from $\tilde T(\tilde x,\;\cdot\;)$, for a given $\tilde x\in\tilde E$, with the choice $(0)$ of $(w_i)_{i\in I}$ above.
Note that an estimator for $\sigma^2(\iota g)$ itself is given by $$\frac{n\sum_{i=1}^n\left|\hat\rho(\hat X_{i-1},\hat Y_i)\right|^2\left|(\iota g)(Y_i)-E_n\right|^2}{\left|\sum_{i=1}^n\hat\rho(\hat X_{i-1},\hat Y_i)\right|^2}\tag8,$$ where $(E_n)_{n\in\mathbb N}$ is an arbitrary estimator of $\hat\mu(\iota g)=\lambda g$) (for example $E_n=I_n\iota g$, n$\in\mathbb N$).
Moreover, it is interesting to note that for a fixed $\tilde x=(i,x')\in\tilde E$, $x:=\varphi_i(x')$ and $\hat x:=(i,x)\in\hat E$, $a_i(x)$ is precisely the variance of the importance sampling estimator of $\lambda f$ using the importance distribution density (with respect to $\lambda$) $$r_{\tilde x}(y):=q_i(y)u'\left(\varphi_i^{-1}(x),\varphi_i^{-1}(y)\right)\;\;\;\text{for }y\in E.$$ So, if $(E_n)_{n\in\mathbb N}$ is any real-valued process with $E_n\xrightarrow{n\to\infty}\lambda f$ almost surely and $(\Sigma_n)_{n\in\mathbb N}$ is an $E$-valued independent process with $\Sigma_n\sim r_{\tilde x}\lambda$ for all $n\in\mathbb N$, then $$\frac1N\sum_{n=1}^N\left|\frac f{r_{\tilde x}}(\Sigma_n)-E_N\right|^2\xrightarrow{N\to\infty}\lambda\frac{|f|^2}{r_{\tilde x}}-|\lambda f|^2=a_i(x)\;\;\;\text{almost surely}\tag{11}.$$ However, from a computational point-of-view, it would be way too cost-intensive to estimate $a_i(x)$ in this way.
Definitions: Let
- $(E,\mathcal E,\lambda)$ be a measure space
- $I$ be a finite nonempty set
- $p,q_i$ be probability densities on $(E,\mathcal E,\lambda)$ for $i\in I$
- $\mu:=p\lambda$
- $w_i:E\to[0,1]$ be $\mathcal E$-measurable for $i\in I$
Assume $$\{p=0\}=\{q_i=0\}\subseteq\{w_i=0\}\;\;\;\text{for all }i\in I\tag1$$ and $$\{p>0\}\subseteq\left\{\sum_{i\in I}w_i=1\right\}\tag2.$$ Now, let
- $(E',\mathcal E',\lambda')$ be a measure space
- $\varphi_i:E'\to E$ be bijective and $(E',\mathcal E)$-measurable with $$\lambda'\circ\varphi_i^{-1}=q_i\lambda\tag3$$ for $i\in I$
- $\zeta$ denote the counting measure on $(I,2^I)$
- $$\tilde t((i,x'),j):=\left.\begin{cases}\left(w_j\circ\varphi_i\right)(x')&\text{, if }(p\circ\varphi_i)(x')>0\\|I|^{-1}&\text{, otherwise}\end{cases}\right\}$$ for $(i,x')\in\tilde E:=I\times E'$ and $j\in I$ and $$\tilde T(\tilde x,\;\cdot\;):=\tilde t(\tilde x,\;\cdot\;)\zeta\;\;\;\text{for }\tilde x\in\tilde E$$
- $u':E'\times E'\to(0,\infty)$ be $\mathcal E'\otimes\mathcal E'$-measurable with $$\int\lambda'({\rm d}y')u'(x',y')=1\;\;\;\text{for all }x'\in E'\tag4$$ and $$U'(x',\;\cdot\;):=u'(x',\;\cdot\;)\lambda'\;\;\;\text{for }x'\in E'$$
- $\hat Q$ denote the transition kernel with density $$\hat q((i,x),(j,y)):=q_j(y)\tilde t\left(\left(i,\varphi_i^{-1}(x)\right),j\right)u'\left(\varphi_j^{-1}(x),\varphi_j^{-1}(y)\right)$$ for $(i,x),(j,y)\in\hat E$ with respect to $\hat\lambda:=\zeta\otimes\lambda$
- $\hat\mu$ denote the measure with density $$\hat p(i,x):=w_i(x)p(x)\;\;\;\text{for }(i,x)\in\hat E$$ with respect to $\hat\lambda$
EDIT:
In my application $E=M^{\{0,\:\ldots\:,\:k\}}$, where $M$ is an oriented 2-dimensional embedded submanifold of $\mathbb R^3$ with surface measure $\sigma_M$ and $k\in\mathbb N$ (in practice, $k\in\{3,\ldots,10\}$). $\mathcal E=\mathcal B(M)^{\otimes\{0,\:\ldots\:,\:k\}}$ and $\lambda=\sigma_M^{\otimes\{0,\:\ldots\:,\:k\}}$. Now $I=\{(s,t):s,t\ge2\text{ and }k+1=s+t\}$ and $$q_{s,\:t}(x)=\rho_1(x_0,x_1)\prod_{i=2}^{t-1}\rho((x_{i-2},x_{i-1}),x_i)\prod_{i=t+2}^k\rho((x_i,x_{i-1}),x_{i-2})\rho_2(x_{k-1},x_k)$$ for all $x\in E$ and $(s,t)\in I$, where $\rho_1,\rho_2$ are probability densities on $\left(M^2,\mathcal B(M)^{\otimes2},\sigma_M^{\otimes2}\right)$ and $\rho:M^2\times M\to[0,\infty)$ is Borel measurable with $$\int\sigma_M({\rm d}z)\rho((x,y),z)=1\;\;\;\text{for all }(x,y)\in M^2.\tag9$$ Moreover, $$p(x)=0.212671f_1(x)+0.715160f_2(x)+0.072169f_3(x)\tag{10}$$ for some $\mathcal E$-measurable $f:E\to[0,\infty)^3$ with $\lambda|f|^2<\infty$ ($\left|\;\cdot\;\right|$ denoting the Euclidean norm). The $g$ in the question is actually $|f|^2$.
Lastly, $(E',\mathcal E',\lambda')=\left([0,1)^d,\mathcal B([0,1))^{\otimes d},\mathcal U_{[0,1)}^{\otimes d}\right)$, $\mathcal U_{[0,1)}$ denoting the uniform distribution on $[0,1)$, where $d:=2k$ and $$u'(x',y'):=\gamma+(1-\gamma)\prod_{i=1}^d\psi_{0,\:\sigma^2}(y'_i-x'_i)\;\;\;\text{for }x',y'\in[0,1)^d,$$ where $\gamma\in[0,1]$, $\sigma>0$ and $\psi_{0,\:\sigma^2}$ is the density of the wrapped normal distribution with mean $0$ and variance $\sigma^2$ on $[0,1]$, i.e. $$\psi_{0,\:\sigma^2}(a)=\sum_{k\in\mathbb Z}\varphi_{0,\:\sigma^2}(k+a)\;\;\;\text{for }a\in[0,1),$$ where $\varphi_{0,\:\sigma^2}$ denotes the density of the normal distribution with mean $0$ and variance $\sigma^2$.
