"Precise" random variables from a distribution

Dear community,

I have an interval from 0 to 1 with a step 0.01 (also can be considered as from 0 to 100). This 101 values unifromly distruted on this interval.

I am trying to randomly draw two values from this interval and expect them to look like 0.15 and 0.25 or 0.9 and 0.56, etc.

I use the following code:

b = 1;
N = 0:0.01:1;
for i = 0:0.01:1
    p1 = rand;
    p2 = rand;
end

This gives me two random values between 0 and 1 but they have a smaller step, for instance, it returns p1 = 0.0538 and p2 = 0.7781.

Any help will be highly appreciated!

Answer 1

[Steven Lord]

Since 0.01 cannot be exactly represented in double precision, your values may not be exactly multiples of one-hundredth. But if you want to draw elements from a specific set of values, you can use randperm (without replacement) or randi (with replacement) to generate indices into that set.

S = [0, 5, 42, -999];
withoutReplacement = randperm(numel(S), 3) % Ask for 3 of the elements of S

withoutReplacement = 1×3

3 2 1

withReplacement = randi(numel(S), [1 3])

withReplacement = 1×3

3 1 3

S(withoutReplacement)

ans = 1×3

42 5 0

S(withReplacement)

ans = 1×3

42 0 42

[Yuriy]

Looks like it generates what I expect, I modified it as follows:

S = 0:0.01:1;
withoutReplacement = randperm(numel(S), 2) % Ask for 3 of the elements of S
withReplacement = randi(numel(S), [1 2])
S(withoutReplacement)
S(withReplacement)

Which produces the following rersut:

withoutReplacement = 92 31

withReplacement = 33 35

ans = 0.9100 0.3000

ans = 0.3200 0.3400

I don't really understand yet what's going on there. Is it possible to express this result ans = 0.3200 0.3400 as two variables like this:

p1 = 0.3200

p2 = 0.3400

and then to use it in my following equations?

[Torsten]

p1 = S(withoutReplacement(1))
p2 = S(withoutReplacement(2))

(same for withReplacement)

[Yuriy]

Great, thanks! Could you please explain me this code? I have a result (please see the screenshot) but don't understand the logic of calculation of S. Why it generates 92, 14? but then it makes it 0.9100 and 0.1300?

S = 0:0.01:1;
withoutReplacement = randperm(numel(S), 2) % Ask for 3 of the elements of S
S(withoutReplacement)
p1 = S(withoutReplacement(1))
p2 = S(withoutReplacement(2))

[Torsten]

Because your S starts at 0, not at 0.01.

Thus element 92 of your array S is 0.91, element 14 of your array S is 0.13.

[Steven Lord]

withoutReplacement is a set of indices that we can use to extract elements of S. It is not the elements of S themselves.

S = [0, 5, 42, -999];
for index = 1:numel(S)
    fprintf("Element %d of S is %d.\n", index, S(index))
end

Element 1 of S is 0.
Element 2 of S is 5.
Element 3 of S is 42.
Element 4 of S is -999.

The first number printed on each of those lines are the indices that would be stored in withoutReplacement. The second number printed on each of those lines are the elements of S stored at those locations.

[Yuriy]

Thank you very much for the explanation! Is there any way to loop it and get an array of all possible combinations and then plot it?

[Torsten]

What do you want to plot if you have a matrix with 101*100 rows and 2 columns for all possible combinations for p1 and p2 ?

This would be the result in your case when taking 2 elements out of 101 (S has 101 elements) without replacement.

[Yuriy]

I want to plot one profit function Profit_1 = p1+p2 for all combinations when p1 < p2. And another profit function Profit_2 = p1 + b, with all possible values of p1 and b ==1. And then see which curve is dominating. I assume that it should look something like on attached drawing.

[Torsten]

Totally unclear how you can arrive at such a plot from the (p1,p2)-matrix. Give me a hint.

What is p ? What is the profit, given p ?

[Yuriy]

Yes, you are right. My apologies. So, I have 2 firms, they draw prices p1 and p2. If p1 < p2, then this firm has two options for a profit:

1. Profit_1 = (p1 + p2)*((1-a)/2 +a) or
2. Profit_2 = (p1+b)(1-a)

Where "a" is just a parameter, any number from interval (0,1), I can set it manually. Equalising it, I can see where they cross each other, but the patterns of these curves aren't clear.

Using formula above for random value I can draw one case and plot it. But I want to plot all possible combinatitons and see a "beahviour" of this curves. This will help me to decide which of two options of profit I have to choose when I set a value for "a".

[Yuriy]

so far I could do like this, but I am not sure it is correct.

a = 0.2;
v = 1;
p_bar = v;
p_low = 0;
p_hat = 0.3;
pd1 = makedist('Uniform','lower',p_low,'upper',p_bar);
pd2 = makedist('Uniform','lower',p_low,'upper',p_bar);


p_1 = (p_low:0.01:p_bar);
p_2 = (p_low:0.01:p_bar);


pdf1 = pdf(pd1,p_1);
cdf1 = cdf(pd1,p_1);

pdf2 = pdf(pd2,p_2);
cdf2 = cdf(pd2,p_2);


f_1 = @(p_1) p_1 * (a + 1-a/2) + p_2 * (a + 1-a/2);
y_1=f_1(p_1);

f_2 = @(p_1) p_1 * (a + 1-a/2) + p_bar*(1-a)/2;
y_2=f_2(p_1);


plot(p_1,y_1,p_2,y_2)