Vectorization of while loop

Hi all,

I have a small function, which generates random numbers in a specific interval. The function tests the generated values in respect to a certain bound and recalculates the random numbers again, if the bound is violated. I have to call this function several times (in this case 60 000 times), so this takes a lot of time (3/4 of the total time). I have read about vectorizing the code to improve performance. But I am totally helpless with this task. Your help would be greatly appreciated.

However, if someone has an idea to rewrite the while loop in a different way that would be great, as well.

Michi

Code:

meanValue = 10;
devValue = 1;
numValues = 10000;
presumption = 1;

for i = 1:60000
    arrayValues = RandomNums(meanValue, devValue, numValues, presumption);
end


function arrayValues = RandomNums (meanValue, devValue, numValues, presumption)
%% create rand nums
% normally distributed with meanValue +/- devValue
    numValues = round(numValues);
    arrayValues = devValue*randn(numValues,1) + meanValue; 
% check confidence-limes for each value and create new one if outside 
    for idx = 1:numel(arrayValues)
        while (arrayValues(idx) > (meanValue + presumption*devValue) || arrayValues(idx) < (meanValue - presumption*devValue))
            arrayValues(idx) = devValue*randn(1,1,'double') + meanValue;
        end
    end
end

Performance:

[Bruno Luong]

You should look for "truncated gaussian" distribution, and posts how to generate them.

[elchico]

Hi Bruno,

thank you for your reply. I've tried that, but as far as I understand, it does not make it better. Maybe, I have missed /missunderstood something?

Thanks again!

Code with Comparison:

meanValue = 10;
devValue = 1;
numValues = 10000;
presumption = 1;

for i = 1:10
    arrayValues = RandomNums(meanValue, devValue, numValues, presumption);
    a = ["arrayValues1",num2str(i)];
    disp(a)
    pause(0.00001);
end
for ii = 1:10
    arrayValues2 = RandomNums2(meanValue, devValue, numValues, presumption);
    a = ["arrayValues2",num2str(ii)];
    disp(a)
    pause(0.00001);
end


function arrayValues = RandomNums (meanValue, devValue, numValues, presumption)
%% create rand nums
% normally distributed with meanValue +/- devValue
    numValues = round(numValues);
    arrayValues = devValue*randn(numValues,1) + meanValue; 
% check confidence-limes for each value and create new one if outside 
    for idx = 1:numel(arrayValues)
        while (arrayValues(idx) > (meanValue + presumption*devValue) || arrayValues(idx) < (meanValue - presumption*devValue))
            arrayValues(idx) = devValue*randn(1,1,'double') + meanValue;
        end
    end
end

function arrayValues2 = RandomNums2 (meanValue, devValue, numValues, presumption)
%% create rand nums
% normally distributed with meanValue +/- devValue
    numValues = round(numValues);
    pd = makedist('Normal','mu',meanValue,'sigma',devValue);
    confidence = presumption*devValue;
    t = truncate(pd,meanValue - confidence,meanValue + confidence);    
    arrayValues2 = random(t,numValues);
end

Performance:

Answer 1

[Bruno Luong]

meanValue = 10;
devValue = 1;
numValues = 1000000;
presumption = 2;
  
tic
% Function from here https://www.mathworks.com/matlabcentral/fileexchange/23832-truncated-gaussian
arrayValues = meanValue + TruncatedGaussian(-devValue, presumption*[-1 1], [1 numValues]);
toc % Elapsed time is 0.055511 seconds for one billions random numbers.

% Check histogram
hist(arrayValues,100)

Histogram obtained

[elchico]

wow, this is so fast.

Thank you!

Can you give me a hint, why it is that fast compared to my initial approach?

[Bruno Luong]

It mainly because it uses different method (non rejection) and inverse error function.

[elchico]

Ok, thank you

Answer 2

[elchico]

Hi Bruno,

one more question to your code: Do you have something similar with Poisson Distribution etc.?

Thanks.

[Bruno Luong]

Unfortunately no.

[elchico]

okay, that is sad for me but: thanks anyways ;-)

May I ask you what is the difference in performance between your code (so fast!) compared to my initial Gauss code? I have to explain this and I am not totally sure about it ...