0
$\begingroup$

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...

$\endgroup$
0
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Yes, I know about erfinv but I don't know how to scale the argument. $\endgroup$ – Pascal Jul 30 '13 at 15:13
  • $\begingroup$ In general, if $\zeta \sim \mathcal{N}(0,1)$, then $z = \mu+\sigma\zeta \sim \mathcal{N}(\mu,\sigma^2)$. $\endgroup$ – Emily Jul 30 '13 at 15:19
  • $\begingroup$ I mean, I'm very confused at what you're asking. At the end of your post, you have an equation for the inverse CDF and you ask "how do I implement this?" You also already apparently know how to do basic change of variables to transform random variables. So, I'm not sure what the "this" is in "how to impement this". $\endgroup$ – Emily Jul 30 '13 at 15:25
0
$\begingroup$

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?

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.