2

we have data that follows he lognormal distribution, knowing it's mean = 5.0163 and standard deviation = 1.0571 we want to generate and draw n (6000) sample within range (3:7.9) that follows the same distribution with monte-carlo method we have this code but it doesn't output samples in the desired range (all samples are smaller than the lower limit)

Data = [ 3 3 3 3 3.3 3.3 3.6 3.9 3.9 3.9 3.9 3.9 3.9 3.9 3.9 3.9 3.9 3.9 3.9 4.2 4.2 4.2 4.2 4.2 4.2 4.5 4.5 4.5 4.5 4.5 4.5 4.5 4.8 4.8 4.8 4.8 4.8 4.8 4.8 4.8 4.8 5.1 5.1 5.1 5.1 5.1 5.1 5.1 5.1 5.1 5.1 5.1 5.1 5.4 5.4 5.4 5.4 5.4 5.4 5.4 5.4 5.4 5.7 6 6 6 6 6 6 6 6 6.3 6.3 6.3 6.3 6.3 6.3 6.3 6.6 6.6 6.6 6.6 6.9 6.9 6.9 7.8];
%draw Lognormal Fit of the data
histfit(Data,[],'lognormal');
% lognormal curve parameters
LN = lognfit(Data);
m=LN(1);
s=LN(2);
% mean and stanard diviation of the associated normal dist.
mu=log(m^2/sqrt(s^2+m^2));
sigma=sqrt(log((s^2/m^2)+1));
%generate random numbers
for i = 1 : 100
X = lognrnd(mu,sigma)
if ((X>=3)&&(X<=7.9))
 X;
end
end
  • 2
    That `if` statement of yours isn't doing anything, whether the statement is true or false. You are calling `X;` which does absolutely naught. – Adriaan Feb 01 '16 at 11:31
  • You can omit that loop by using `X = lognrnd(mu,sigma,100,1)` if you want to have 100 values. Btw: mu = 0.45, sigma = 0.136 in your example code. – Matthias W. Feb 01 '16 at 11:34
  • [Take a look at the documentation for logical indexing, first example](http://de.mathworks.com/company/newsletters/articles/matrix-indexing-in-matlab.html?refresh=true) – Daniel Feb 01 '16 at 12:10
  • @Matthias, you are right mu and sigma weren't correct. now we can generate the 6000 samples using 'X = lognrnd(mu,sigma,6000,1)'. What about limiting the generated values to b in the range from 3 to 7.9 – user3332603 Feb 01 '16 at 12:36
  • @user3332603 you can "filter" your results using logical indexing: `X( (X >= 3) & (X <= 7.9) )`, but this will obviously drop the values that are not within your given interval - giving you less than the 100 values. – Matthias W. Feb 01 '16 at 12:45
  • 1
    Possible duplicates or references: [SO - How to generate random numbers of lognormal distribution within specific range in Matlab](http://stackoverflow.com/questions/3070183/how-to-generate-random-numbers-of-lognormal-distribution-within-specific-range-i), [MATLAB Central - Log normal distribution between two value](http://www.mathworks.com/matlabcentral/newsreader/view_thread/337511) – Matthias W. Feb 01 '16 at 12:51

1 Answers1

0

It is not possible to guarantee samples from a Lognormal distribution will fall within a certain (finite) range as the support for the distribution is unbounded (all positive real numbers). However, we can invoke MATLAB's truncate function to easily deal with this issue (see below).

Parameterization: Note that if X ~ Normal(mu,sigma), then if Y = ln(X), Y ~ Lognormal(mu,sigma). Thus, the Lognormal is often parameterized by the associated Normal distribution. The mean, E[Y] is not equal to mu.

This is illustrated by the code below.

mu = 5; 
sigma = 1;
pd = makedist('Lognormal',mu,sigma)
Y = random(pd,50000,1);                 % 50000 Lognormal samples
X = log(Y);  

>> mean(X)                              
ans =
    4.9922
>> std(X)
ans =
    1.0007

Restriction to the range [a,b]:
Given the parameters mu and sigma obtained from fitting the Lognormal, then restriction to an given interval, [a,b] is simple with MATLAB's probability distribution objects (introduced in R2013a).

mu = 2.5;       % Assume mu & sigma found from fitting 
sigma = 1.4;
pd = makedist('Lognormal',mu,sigma);    % Creates the Lognormal distribution object

a = 10;
b = 100;
N = 10000;
pdT = truncate(pd,a,b);              % Truncates to the interval [a,b]

Y = random(pdT,N,1);                 % Generate samples from truncated Lognormal

The difference is clear.

Lognormal & Truncated Lognormal

W = random(pd,N,1);                  % Generate samples from Lognormal 

figure
s(1) = subplot(2,1,1)
histogram(Y,'Normalization','pdf')
xlim([0,125])
s(2) = subplot(2,1,2)
histogram(W,'Normalization','pdf','FaceColor','g','FaceAlpha',0.35)
xlim([0 125])
title(s(2),'Lognormal(\mu=2.5,\sigma=1.4)')
title(s(1),'Lognormal(\mu=2.5,\sigma=1.4) truncated to [10,100]')
SecretAgentMan
  • 2,856
  • 7
  • 21
  • 41