Quasi Monte Carlo estimation of logit-normal density integrals I am considering the integral
$$ I(y \mid \mu, \sigma) = \int_y^1  \frac{\exp \left\{ \frac{-1}{2\sigma^2}[\textrm{logit}(x)-\mu)]^2 \right\}}{\sigma \sqrt{2\pi} (1-x)}\textrm{d}x$$
which for $y=0$ is the expectation of a random variable $X$ such that $X \sim \textrm{logit-Normal}(\mu, \sigma^2)$. Unfortunately, the integral is intractable so numerical methods are required. This integral is a step I need to take as part of an optimisation routine so I require a deterministic and fast estimate of $I(y)$. For unfortunate choices of $\mu$, $\sigma$ this integral can be difficult to calculate numerically since the density imight infinite at $0$ or $1$.
I know that a quasi Monte Carlo (QMC) estimate of $I(0)$ is
$$I(0) = E(X) \approx \frac{1}{K-1}\sum_{i=1}^{K-1} G \left(\Phi^{-1}_{\mu, \sigma^2}(i/K)\right)$$
with $G(x) = 1/(1 + e^{-x})$ and $\Phi^{-1}_{\mu, \sigma^2}(\cdot)$ is the inverse of a normal cdf with mean $\mu$ and variance $\sigma^2$. This estimate is determinsitic and fast! However I'm struggling to generalise to $I(y)$ with $0<y<1$
I have tried numerical integration in $\verb|R|$ using $\verb|integrate()|$ and for the most part is great, but can be unstable for unfortunate choices of $\mu$, $\sigma$, hence QMC.
I have also tried starting the sum at the first value of $i$ such that $i>Ky$ but this gives a really bad estimate (I adjusted the denominator to account for the 'length' of the sum changing).
For instance, a regular Monte Carlo estimate of $I(0.8;1,1) \approx 0.6397$ whereas using a QMC estimate with the sum starting at the first value of $i$ such that $i>Ky$ is $I(0.8 \mid 1, 1) \approx 0.9106$ which is clearly way off. I verified by MC estimate is valid by comparing to a vanilla numerical integration. Comparing regular MC estimates of $I(0)$ to QMC estimates of $I(0)$ also give good agreement; there's something I don't understand about this 'truncation'.
Failing this, are there alternate, fast methods of computing $I(y|\mu, \sigma)$?
There's some $\verb|R|$ code below which outlines the two computations I illustrated
## regular MC
X <- rnorm(10000,1,1)
X <- mean((expit(X))*(expit(X)>1))
> X
[1] 0.6397533

## QMC
Z <- ((1:(10000-1))/10000)
Z <- Z[Z > 0.8]
mean(expit(qnorm( Z ,1,1)))
[1] 0.9106322



 A: The way I understand the integral, it writes as
\begin{align*}
I(y \mid \mu, \sigma) &= \int_y^1  \frac{\exp \left\{ \frac{-1}{2\sigma^2}[\textrm{logit}(x)-\mu]^2 \right\}}{\sigma \sqrt{2\pi} (1-x)}\textrm{d}x\\
&= \int_y^1  \frac{\varphi_{\mu,\sigma^2}(\textrm{logit}(x))}{(1-x)}\textrm{d}x\\
&\stackrel{z(x)=\textrm{logit}(x)}{=} \int_y^1  \varphi_{\mu,\sigma^2}(z(x)) x\frac{\textrm{d}z(x)}{\textrm{d}x}^\text{(i)}\textrm{d}x\\
&= \int_{\textrm{logit}(y)}^\infty G(z)^\text{(ii)} \varphi_{\mu,\sigma^2}(z)\textrm{d}z\\
\end{align*}
since
(i). the logit function $z(x)=\textrm{logit}(x)=\log\{x/(1-x)\}$ has derivative$$\frac{\textrm{d}z(x)}{\textrm{d}x}=\frac{1}{x(1-x)^2}$$
(ii) the inverse of the logit function is $x(z)=e^x/(1+e^x)=1/(1+e^{-x})=G(z)$
Hence it is a Gaussian expectation and a deterministic approximation of the integral should be
$$\frac{1}{K-1}\sum_{i=1}^{K-1} G \left(\Phi^{-1}_{\mu, \sigma^2}(i/K)\right)
\times\mathbb I_{\underbrace{\textrm{logit}(y)\le\Phi^{-1}_{\mu, \sigma^2}[i/K)]}_{\Phi_{\mu, \sigma^2}[\textrm{logit}(y)]\le i/K}}$$
I would however not call it a quasi-Monte Carlo approximation, but rather a Riemann sum approximation.
