Shift interval of log-normally distributed latin hypercube samples

Question

first of all I'm not sure if this part of StackExchange is the right one because my question is mainly on a way to implement something in MATLAB.

Ok, now let me try to pack my whole question in one sentence:

How do I implement latin hypercube sampling of log-normal distribution whose interval is shifted from $x\in (0,+\infty)$ to $x\in (1,+\infty)$?

More detailed:

I have a fit algorithm with a bunch of variables that are defined on different intervals ($x_i\in[0,1]$, $y_i\in[0,+\infty)$ and $z_i\in[1,+\infty)$). The variables $x_i$ are normally distributed while $y_i$ and $z_i$ should be log-normally distributed. I now want to sample a couple of starting values for the fit algorithm for each of these values from a latin hypercube according to their distribution.

For that I found the MATLAB package latin hypercube sampling whose function latin_hs (link) can be used to tackle the $x_i$s.

For the variables $y_i$ I modified latin_hs in the following way:

function s = lhslogn(xmean, xsd, nsample, nvar)

r = rand(nsample, nvar);
s = zeros(nsample, nvar);

for i = 1:nvar
    idx = randperm(nsample);
    P = (idx' - r(:,i))/nsample;
    s(:,i) = logninv(P, xmean(i), xsd(i));
end

where logninv a built-in function calculating the inverse cummulative distribution function.

Now, my question is: what do I have to do to sample the variables $z_i$?

---

Notes:

If I understood it correctly, I can just use $y = x-1$ to shift the log-normal probability density function to the right by $1$, cf.:

$$ P(y;\mu,\sigma) = \frac{1}{\sqrt{2\pi}\sigma y}\exp\left(-\frac{(\ln(y)-\mu)^2}{2\sigma^2}\right) $$

which then gives the cumulative distribution function

$$ C(y; \mu, \sigma) = \frac{1}{2}\left(1+\frac{\ln(y)-\mu}{\sqrt 2 \sigma}\right) = \Phi\left(\frac{\ln(y)-\mu}{\sigma}\right) $$

where $\Phi$ is the CDF of the standard normal distribution.

The inverse CDF can now be calculated as

$$ D(P) = \exp\left(\mu+\sigma\Phi^{-1}(P)\right). $$

But I now have absolutely no idea how to implement this...

Answer 1

Ok, I think I figured it out.

The part of the previous code

s(:,i) = logninv(P, xmean(i), xsd(i));

has to be replaced with

logx0 = -sqrt(2).*erfcinv(2*P) + 1;
s(:,i) = exp(xsd(i).*logx0 + xmean(i));

The important part here is the +1 in the first line.

Explanation: Basically all I want is the corresponding x-values to some specific y-values of the log-normal CDF. This can be done by inverting the CDF. The inverse of the normal CDF is known as the probit function which is

$$ \text{probit}(P) = -\sqrt{2}\text{erfc}^{-1}(2P). $$

A shift of the CDF to the right by 1 implicates a up-shift of the probit function by 1. That's why I simply added the +1 to the probit function.

Is that correct?

Answer 2

The function erfinv will give you the inverse CDF $\Phi^{-1}$. Simply scale your argument appropriately.