Loop to pick random variables from separate probability distributions

I have a dataset that can be described by a beta probability distribution (alpha = 3.18, beta = 44.87). What I want to do is pick a random value from that beta probability distribution, and depending how many standard deviations away from the mean that value is, I would like to select a second random value from a lognormal distrbution (mu=45.89, sigma=151.595) that is also within the specific standard deviation range.

Summarized: if a value of 0.17 is picked from my beta distrbution (which lies to the right of the mean and is greater than 2nd std, but less than the 3rd std) then I want to pick a random value from my lognormal distrbution that is also to the right of the mean and greater than the 2nd std but less than the 3rd std).

Firstly, I apologize as I am very new to this (1 week in). But the main issue with my code is that I do not properly know how to pick a random value from my lognormal distrbution within a specific standard deviation range. In the below code I have this step simply as the placeholder: Perm = lognrnd(45.89,151.5956,[1 1]). I am also missing one line of code at the end of my for loop which would copy my porRandom value into one column and Perm into a second column after each of the 300 runs, resulting in a 300x2 matrix.

Thanks for taking the time!

number_of_runs = 300;                                                                                  % Run this loop 300 times
porRandom = betarnd(3.1800,44.87,[1 1]);                                                               % porRandom = random number generated from a beta distribution within a 1x1 matrix
Perm = (number_of_runs,1);                                                                             % Generate a value for perm for each loop

for n = 1:number_of_runs;                                                                   
if porRandom <= (mean(porDist) + std(porDist)) OR >= (mean(porDist) - std(porDist));                   % if porRandom is within one standard deviation of the mean of a beta distrbution (porDist)
    Perm = lognrnd(45.89,151.5956,[1 1])                                                               % then let Perm be equal to a random value within one standard deviation of a defined lognormal distrbution
else if porRandom < (mean(porDist) - (std(porDist) AND > (mean(porDist) - ((std(porDist))*2));         % else if porRandom is within a 1st to 2nd standard deviation range (less than the mean)
        Perm = lognrnd(45.89,151.5956,[1 1])                                                           % then let Perm be equal to a random value within a 1st to 2nd standard deviation range (less than the mean) of a defined lognormal distrbution
    else if porRandom < (mean(porDist) - (((std(porDist))*2))                                          % else if porRandom is greater than 2 standard deviations (less than the mean)
            Perm = lognrnd(45.89,151.5956,[1 1]                                                        % then let Perm be equal to a random value that is greater than 2 standard deviations (less than the mean)
        else if porRandom > (mean(porDist) + (std(porDist) AND < (mean(porDist) + ((std(porDist))*2)); % else if porRandom is within a 1st to 2nd standard deviation range (greater than the mean)
                Perm = lognrnd(45.89,151.5956,[1 1]                                                    % then let Perm be equal to a random value within a 1st to 2nd standard deviation range (greater than the mean) of a defined lognormal distrbution 
            else porRandom > (mean(porDist) + (((std(porDist))*2));                                    % else porRandom is greater than two standard deviations (greater than the mean)
                    Perm = lognrnd(45.89,151.5956,[1 1]                                                % then let Perm be equal to a random value that is greater than 2 standard deviations (greater than the mean) 
end                                                                                        

Answer 1

[Jeff Miller]

It sounds like you want to select the second number from a truncated version of the lognormal distribution. If I understand what you are after, here is one way to do it with Cupid (just generating one random number for an example):

% Set up the beta distribution
% dist1 = Beta(3.18,44.87);
dist1.PlotDens   % Make sure this is the distribution you want
dist1sd = dist1.SD

% Set up the lognormal
% dist2 = LognormalMS(45.89,151.595);
dist2.PlotDens  % Make sure this is the distribution you want
dist2mn = dist2.Mean
dist2sd = dist2.SD

% Generate a random beta and find out where it lies in terms of SD's
dist1rnd = 0.17;  % or dist1rnd = dist1.Random to pick this randomly
dist1lowsd = floor( (dist1rnd-dist1.Mean) / dist1sd )
dist1hisd = ceil( (dist1rnd-dist1.Mean) / dist1sd )

% Bound the lognormal to lie within the same range of sds.
boundedlognormal = TruncatedX(dist2,dist2mn+dist1lowsd*dist2sd,dist2mn+dist1hisd*dist2sd);
desiredrandom = boundedlognormal.Random   % Get a random lognormal from that restricted range
desiredrandomsd = (desiredrandom - dist2mn) / dist2sd  % Check: Should be between dist1lowsd & dist1hisd

[Evan Renaud]

Hi Jeff,

Thank you for the reply. Yep this is exactly what I want. In the time between now and when I posted this I also figured out that what I needed to look at was truncated distributions. Not that I don't want to use your code, but since I am new at this I think it would be good practice for me to build something too. I have updated my code a bit, and I was curious if you could also take a look at it:

pdPor = beta distribution

pdPerm = lognormal distribution

number_of_runs == 300;                                                                                  % Run this loop 300 times
porRandom == betarnd(3.1800,44.87,[1 1];                                                               % porRandom = random number generated from a beta distribution
Perm == number_of_runs,1;                                                                             % Generate a value for perm for each loop

for n = 1:number_of_runs;                                                                   
if porRandom <= (mean(pdPor) + std(pdPor)) OR >= (mean(pdPor) - std(pdPor));                           % if porRandom is within one standard deviation of the mean
               t = truncate(pdPerm,-1.56163,3.587676117)                                               % then truncate pdPerm distribution to only include data within 1 standard deviation of the mean
               Perm = random(t,1.01302,2.579);                                                         % then let Perm be equal to a random value within the truncated distribution, t
else if porRandom < (mean(pdPor) - (std(pdPor) AND > (mean(pdPor) - ((std(pdPor))*2));                 % else if porRandom is within a 1st to 2nd standard deviation range (less than the mean)
               t = truncate(pdPerm,-4.13628,-1.56163)                                                  % then truncate pdPerm distribution to only include the range of 1 and 2 standard deviations less than the mean
               Perm = random(t,1.01302,2.579);                                                         % then let Perm be equal to a random value within the truncated distribution, t
    else if porRandom < (mean(pdPor) - (((std(pdPor))*2))                                              % else if porRandom is greater than 2 standard deviations (less than the mean)
               t = truncate(pdPerm,inf,-4.13628)                                                       % then truncate pdPerm distribution to only include the range of 2 standard deviations to infinity less than the mean
               Perm = random(t,1.01302,2.579);                                                         % then let Perm be equal to a random value within the truncated distribution, t
        else if porRandom > (mean(pdPor) + (std(pdPor) AND < (mean(pdPor) + ((std(pdPor))*2));         % else if porRandom is within a 1st to 2nd standard deviation range (greater than the mean)
               t = truncate(pdPerm,3.587676,4.574652)                                                  % then truncate pdPerm distribution to only include the range of 1 to 2 standard deviations greater than the mean
               Perm = random(t,1.01302,2.579);                                                         % then let Perm be equal to a random value within the truncated distribution, t 
            else porRandom > (mean(pdPor) + (((std(pdPor))*2));                                        % else porRandom is greater than two standard deviations (greater than the mean)
               t = truncate(pdPerm,4.574652,inf)                                                       % then truncate pdPerm distrbution to only include the range of 2 standard deviations to infinity greater than the mean
               Perm = random(t,1.01302,2.579);                                                         % then let Perm be equal to a random value within the truncated distribution, t
end                                                                                        

[Jeff Miller]

Well, I'm not exactly sure what you are trying to achieve here (partly because there are still a lot of syntax problems), but a few comments...

Notice that porRandom is only generated once, before the loop starts. This means that all 300 of the Perm values will be generated by the same one of the if/elseif clauses (i.e., within the same range of the truncated lognormal). I interpreted the original question to mean that you wanted to generate (porRandom,Perm) pairs where porRandom was truly random and Perm was constrained with your SD idea.

Consider the first part of your if statement, which I would rewrite

if porRandom <= (mean(pdPor) + std(pdPor)) || ...
        porRandom >= (mean(pdPor) - std(pdPor))
end

I think that is always going to be true. Don't you really want && instead of || here?

If speed is an issue, you could improve this by forming the different truncated distributions in advance and then just picking the one you want in each of the 300 iterations.

Sorry, I don't know really know how truncate and random work together, but "random(t,1.01302,2.579)" looks wrong to me. I wouldn't think you would pass the parameters of t into random this way. Instead, I'd think random would take t and then the dimensions of the output array of random numbers to be produced.