2

I have attempted to generate a triangular probability distribution in Matlab, but was not successful. I used the formula at http://en.wikipedia.org/wiki/Triangular_distribution.

n = 10000000;

a = 0.2;
b = 0.7;
c = 0.5;

u = sqrt(rand(n, 1));

x = zeros(n, 1);
for i = 1:n
    U = u(i);
    if U < (c-a)/(b-a)
        X = a + sqrt(U*(b-a)*(c-a));
    else
        X = b - sqrt((1-U)*(b-a)*(b-c));        
    end
    x(i) = X;
end

hist(x, 100);

The histogram looks like so:

enter image description here

Doesn't look like much of a triangle to me. What's the problem? Am I abusing rand(n)?

Superbest
  • 25,318
  • 14
  • 62
  • 134

3 Answers3

5

you can add up two uniform distributions, the distribution graphs convolve, and you get a triangular distribution.

easy-to-understand example: rolling two dice, each action has uniform distribution to result in a number from 1-6, combined action has triangular distribution to result in a number 2-12

edit: minimal working example:

a=randint(10000,1,10);
b=randint(10000,1,10);

c=a+b;

hist(c,max(c)-min(c)+1)

edit2: looked in your script again. It's working but you've made one mistake:

u = sqrt(rand(n, 1));

should be

u = rand(n, 1);

edit3: optimized code

n = 10000000;

a = 0.2;
b = 0.7;
c = 0.5;

u = rand(n, 1);
x = zeros(n, 1);

idx = find(u < (c-a)/(b-a));
x(idx) = a + sqrt(u(idx)*(b-a)*(c-a));
idx =setdiff(1:n,idx);
x(idx) = b - sqrt((1-u(idx))*(b-a)*(b-c));
hist(x, 100);

Gunther Struyf
  • 11,158
  • 2
  • 34
  • 58
  • 1
    The sum of two uniforms is by far the easiest way to do it. However, you can also invert the CDF, which requires only one rand call, but slightly more math. –  Feb 11 '12 at 17:41
  • 2
    yeah, the sum of two uniforms is easiest, but it is only usefull for a symmetrical triangular distribution. The code in topicstart allows for any triangular shape. btw I believe you can still optimize that code by using logical indexing to remove the loop.. – Gunther Struyf Feb 11 '12 at 19:29
  • 2
    I am not sure why you would call "edit3: optimized code" optimized as it takes twice as much time run than the OP's (corrected) code. Vectorizing a simple for-loop using 'find' and 'setdiff' does not make the code optimal. – Kavka Feb 12 '12 at 16:28
3

This example uses the makedist() and pdf() commands.

a = 2; m = 7; b = 10;
N = 50000;                             % Number of samples
pd = makedist('Triangular',a,m,b);     % Create probability distribution object

T = random(pd,N,1);                    % Generate samples from distribution

Triangular Distribution with lowerbound a = 7, mode m = 10, and upperbound b = 10.
Triangular Distribution

% Plot PDF & Compare with Generated Sample
X = (a-2:.1:b+2);

figure, hold on, box on
histogram(T,'Normalization','pdf') % Note normalization-pdf option name-value pair
title([num2str(N) ' Samples'])
plot(X,pdf(pd,X),'r--','LineWidth',1.8)
legend('Empirical Density','Theoretical Density','Location','northwest')

MATLAB introduced makedist() in R2013a. Requires Stats toolbox.

Reference:
Triangular Distribution

SecretAgentMan
  • 2,856
  • 7
  • 21
  • 41
1

Change

u = sqrt(rand(n, 1));

to

u = rand(n, 1);

The nice thing about this formula is that you can distribute a sample from a general triangle distribution with a single random sample.

Peter de Rivaz
  • 33,126
  • 4
  • 46
  • 75