lognrnd Function does not work properly at high variance

I am trying to generate i.i.d lognormal random variables X_i, under the constraint of the average power is one, i.e., E(X^2)=1. However, when using this "lognrnd" function and check the average power of large number of realizations, the output does comply to the constraint of power. I noticed that the function works well at small values of variance. However when increasing the variance the function does not work properly.

   v = .9999; % in order to get sigma=3
m = sqrt(1-v); % normalizing the power of the generated lognormal samples 
mu = log((m^2)/sqrt(v+m^2));
sigma = sqrt(log(v/(m^2)+1))
[M,V]= lognstat(mu,sigma);
X = lognrnd(mu,sigma,1,1e6);
MX = mean(X)
VX = var(X)

the answer

sigma =

    3.0349

MX =

    0.0103

VX =

    0.1825

the sum of VX, MX^2 is not equal to one. !!!!

Answer 1

[the cyclist]

I suspect there is no problem at high variance -- but that you would need an astronomically large number of samples to verify the empirical power is equal to 1 (on average).

Here is a less extreme example. I ran your code, except that I set v = 0.9. I actually ran it 5000 times, and stored the resulting value of P = VX + MX^2 from each run. Each time, P is a little different from 1 -- sometimes quite different. The mean value of P was about 0.999, and the variance of P was about 0.005, so the values are fairly tightly clustered around 1. Here's what the distribution of P looked like:

Now look what happens when I set v = 0.99.

Notice a few things. The extent of the x-axis is much greater, because that's needed to capture the whole distribution (which now has variance 5.2). I am much more likely to get a value far from 1 because of the wide distribution. Also, because the distribution is also skewed, I am much more likely to get a value below 1 -- but I might get one much greater than 1. (However, the mean value is about 0.993, very close to 1.)

This effect will get more and more exaggerated as v gets large, leading to the result you found.

[Ahmed Shaban]

Thanks so much for you effort and help. I checked the two figures and I think that You concluded the same result I had before. However, I got a big difference between the variance of the generated data-even averaging over large number 10e6- and the meant variance.

Do you suggest any modification in order to get the good result?

[the cyclist]

Frankly, I'm not sure what you expect, or what you mean by a "good result". The facts are the facts.

Here is another run, where I used v=0.9999, increased the number of samples from 1e6 to 1e8, and calculated the empirical power 500 times.

You see that even with 100 times larger sampling, the mode is pretty far from 1. These are immutable facts. 1e6 is just not a "large number" of iterations for your distribution.

Answer 2

[Roger Stafford]

If you look at the probability density function that corresponds to your particular mu and sigma combination, you will see that a huge number of samples would have to be placed extremely close to zero by ‘lognrnd’ in order to achieve any reasonable approximation to your requested lognormal distribution. This in turn means that instead of doing only a million samples you would probably need to do billions - maybe even billions of billions - of samples to achieve a good approximation. In other words, you have given ‘lognrnd’ a most awkward set of parameters and it requires a great many more than a mere million samples to achieve the results of the sample variance you are trying for. Look at it from the point of view of your m and v. The mean should be a tiny .01, all samples must be positive, and yet you are requiring a variance of almost 1. That is difficult to achieve for almost any kind of distribution.

The formula for the pdf is:

   f(x) = 1/(x*sigma*sqrt(2*pi)).*exp(-(log(x)-mu).^2/(2*sigma^2));