Calculating and plotting conditional distribution.

Dear community,

I hope you are all good here. Could you please help me with a code for the next problem. I have the following code

lambda = 0.3;
p_bar = 1;
p_tilda = (1 - lambda)/(1 + lambda);
S = 0:0.01:p_tilda;
withoutReplacement = randperm(numel(S), 2); % Ask for 2 of the elements of S
S(withoutReplacement);
p1 = S(withoutReplacement(1));
p2 = S(withoutReplacement(2));

if p1 < p2 
pe1 = @(p1,p2) p1 * ((1 - lambda)/2 + lambda)  + p2 * ((1 - lambda)/2 + lambda);
pe2 = @(p1) p1 * ((1 - lambda)/2 + lambda) + p_bar * (1 - lambda)/2;
y1 = pe1(p1,p2);
y2 = pe2(p1);
else 
    y1 = 0;
    y2 = 0;
end

I generate two random values from the interval S and then calculate two functions if p1 < p2.

Now I want to calculate and plot a distribution for all possible values of y1 and y2 when p1 < p2. I assume that I should make a loop of many drawings for all p1 < p2, but my MathLab skills is not sufficient to do it. I hope someone knows how to code it. Cheers!

Answer 1

[Torsten]

You mean this

lambda = 0.3;
p_bar = 1;
p_tilda = (1 - lambda)/(1 + lambda);
S = 0:0.01:p_tilda;
n = 1000000;
icount = 0;
for i = 1:n
    withoutReplacement = randperm(numel(S), 2); % Ask for 3 of the elements of S
    p1 = S(withoutReplacement(1));
    p2 = S(withoutReplacement(2));
    if p1 < p2 
        icount = icount + 1;
        y1(icount) = p1 * ((1 - lambda)/2 + lambda)  + p2 * ((1 - lambda)/2 + lambda);
        y2(icount) = p1 * ((1 - lambda)/2 + lambda) + p_bar * (1 - lambda)/2;
        %else 
        %    y1(i) = 0;
        %    y2(i) = 0;
    end 
end
figure(1)
ksdensity(y1)
hold on
ksdensity(y2)
hold off

or this

lambda = 0.3;
p_bar = 1;
p_tilda = (1 - lambda)/(1 + lambda);
S = 0:0.01:p_tilda;
n = 10000000;
for i = 1:n
    withoutReplacement = randperm(numel(S), 2); % Ask for 3 of the elements of S
    p1 = S(withoutReplacement(1));
    p2 = S(withoutReplacement(2));
    if p1 < p2 
        y1(i) = p1 * ((1 - lambda)/2 + lambda)  + p2 * ((1 - lambda)/2 + lambda);
        y2(i) = p1 * ((1 - lambda)/2 + lambda) + p_bar * (1 - lambda)/2;
    else 
        y1(i) = 0;
        y2(i) = 0;
    end 
end
figure(2)
ksdensity(y1)
hold on
ksdensity(y2)
hold off

?

[Yuriy]

Hi Torsten! Thank you for the suggestions. I don't think it correctly illustrates what I am trying to plot. From my code you can see that for all p1 < p2 profits y1 < y2. So, I think that the plot should look like two inverted U-shape distributions where one always dominates another for all possible combinations when p1 < p2.

I tried to slighlty modify your code like this:

lambda = 0.3;
p_bar = 1;
p_tilda = (1 - lambda)/(1 + lambda);
S = 0:0.01:p_tilda;
n = 1000000;
for i = 1:n
    withoutReplacement = randperm(numel(S), 2); % Ask for 2 of the elements of S
    p1 = S(withoutReplacement(1));
    p2 = S(withoutReplacement(2));
    if p1 < p2 
        y1 = p1 * ((1 - lambda)/2 + lambda)  + p2 * ((1 - lambda)/2 + lambda);
        y2 = p1 * ((1 - lambda)/2 + lambda) + p_bar * (1 - lambda)/2;
        else 
            y1 = 0;
            y2 = 0;
    end 
end
figure(1)
ksdensity(y1)
hold on
ksdensity(y2)
hold off

And this looks more or less closer to my expecteations but I am not sure why it has negative values, because p1, p2, y1, y2 values are strrictly positive. I'd expect y1 curve be strictly dominating on the interval (0, p_tilda). Maybe you can tell what I am doing wrong here. Thanks a lot!

[Torsten]

I don't think it correctly illustrates what I am trying to plot. From my code you can see that for all p1 < p2 profits y1 < y2.

The plot of the two curves is an estimate for the probability that y1 / y2 attain a certain value. The x-axis is this value, the y axis is in principle the probability that this value results from the random experiment. It has nothing to do with a plot of y over p.

And this looks more or less closer to my expecteations but I am not sure why it has negative values, because p1, p2, y1, y2 values are strrictly positive.

After changing the for-loop as you did, y1 and y2 are no longer arrays, but you overwrite the values they attain in the loop every time so that in the end, y1 and y2 are two single values. And from these two single values you try to estimate a probability density curve. Total nonsense.

[Yuriy]

Thank you very much for the explanation! Could you tell me please is it possible to plot y over p having y1 and y2 as arrays like in the code which you suggest above?

[Torsten]

That's where we were already: give me a hint how you imagine to plot a matrix of the form (p1,p2,y1) with n = 1000000 rows and p1, p2 randomly chosen.

Plotting y2 over p1 is of course easy:

lambda = 0.3;
p_bar = 1;
p_tilda = (1 - lambda)/(1 + lambda);
S = 0:0.01:p_tilda;
y2 = S * ((1 - lambda)/2 + lambda) + p_bar * (1 - lambda)/2;
plot(S,y2)

For y1 you can get a similar plot, but over the triangle (p1,p2), 0 <= p1 < p2 <= 1, as a two-dimensional surface plot.

[Torsten]

Or if you use the smallest and highest possible value for p2, given p1, for y1, the plots look like:

lambda = 0.3;
p_bar = 1;
p_tilda = (1 - lambda)/(1 + lambda);
S = 0:0.01:p_tilda;
y1_low = 2*S * ((1 - lambda)/2 + lambda) ;
y1_high = (S + p_tilda) * ((1 - lambda)/2 + lambda);
plot(S,[y1_low;y1_high])

Note that y1_high = y2.

[Yuriy]

Thanks a lot! I think the best I can do is to generate all possible pairs of values for p1<p2 in a separate file, like set of pairs in the rows and values of p1 and p2 in two columns. And then use this file to calculate all possible values of y1 and y2. And then to plot it. Now I will try to find how to generate such a table in MatLab, if it's possible, of course.

Thanks a lot!

[Torsten]

As said, the problem is not to generate the table, but to find a way to plot it sensefully.

[Yuriy]

Thank you! I understood. So, I did this file with all possible variations of p1 < p2 and then made a plot in excel (I've attached the file and the graph. Not sure if graph is correct and sensible though, but it looks to me like I managed to plot all possible values of Y against P and can see which is dominating). I wonder is it possible to plot y1 and y2 against p using this file in Matlab?

[Torsten]

As said (in the comment you deleted): I don't understand the Excel graph - so I also cannot help you to plot it in MATLAB.

And didn't you claim that y1 should always be smaller or equal y2 ? In your graph, this doesn't seem to be true.

[Yuriy]

I ploted for the whole interval (p_low, p_bar). On the interval (p_low, p_tilda) Y1 is always lower than Y2 for all p1 < p2. But on the interval (p_tilda, p_bar) it is not the case anymore.

[Torsten]

In your MATLAB description, p1 and p2 can only be drawn from (p_low,p_tilda).

So you changed the rules in Excel ?

[Yuriy]

If I will keep only values on the interval (p_low, p_tilda) then the plot will look like this:

[Torsten]

It you connect the lower start points of the blue lines, you get (together with the red line) exactly the lower and upper curves for y1 I plotted.

[Yuriy]

Ah yes, I can see it now. Thanks a lot! I assume, it means that my plot is correct and makes sense too. However, is it possible to plot it in MatLab? With or without the file I used? And how CDF will look then?

[Torsten]

lambda = 0.3;
p_bar = 1;
p_tilda = (1 - lambda)/(1 + lambda);
S = 0:0.01:p_tilda;
icount = 0;
for i = 1:numel(S)-1
    p1 = S(i);
    for j = i+1:numel(S)
        p2 = S(j);
        icount = icount + 1;
        y1(icount) = p1 * ((1 - lambda)/2 + lambda)  + p2 * ((1 - lambda)/2 + lambda);
        y2(icount) = p1 * ((1 - lambda)/2 + lambda)  + p_bar * (1 - lambda)/2;
    end
end
figure(1)
plot(1:icount,y1)
hold on
plot(1:icount,y2)
hold off

figure(2)
ecdf(y1)
hold on
ecdf(y2)
hold off

[Yuriy]

Thanks a lot! Could you please explain me the code a little? I can see that p1 and p2 aren't the arrays there.

[Torsten]

There is nothing special about the plot. I built all (p1,p2) combinations with p1 < p2, computed y1 and y2 and plotted y1 and y2 over the combination number.

The number of combinations is

sum_{i=1}^{i=100} i = 100*101/2 = 5050.

[Yuriy]

awesome! thanks! I was trying to do the same in Wolfram Alpha, because there you can make a really good 3D live graphs. I wonder is it possible to do the same in MatLab? To plot it for 3 axis p1, p2, and y. Like on attached drawing. Maybe I need to create a separate question for it?

[Torsten]

At the moment, I don't see a simple solution since you don't have to plot over a rectangular domain, but a triangle p1 < p2.

I'd open a new question.

Say in your text that with the code

lambda = 0.3;
p_bar = 1;
p_tilda = (1 - lambda)/(1 + lambda);
S = 0:0.01:p_tilda;
[S1 S2] = ndgrid(S,S);
y1 = S1 * ((1 - lambda)/2 + lambda)  + S2 * ((1 - lambda)/2 + lambda);
y2 = S1 * ((1 - lambda)/2 + lambda)  + p_bar * (1 - lambda)/2;
figure(1)
surf(S1,S2,y1)

figure(2)
surf(S1,S2,y2)

you want to make a surface plot only over S1 < S2.

[Torsten]

lambda = 0.3;
p_bar = 1;
p_tilda = (1 - lambda)/(1 + lambda);
S = 0:0.01:p_tilda;
[S1,S2] = ndgrid(S,S);
idx = find(S1 < S2);
S1 = S1(idx);
S2 = S2(idx);
S = [S1,S2];
T = delaunayTriangulation(S);
y1 = T.Points(:,1) * ((1 - lambda)/2 + lambda)  + T.Points(:,2) * ((1 - lambda)/2 + lambda);
y2 = T.Points(:,1) * ((1 - lambda)/2 + lambda)  + p_bar * (1 - lambda)/2;
figure(1)
trisurf(T.ConnectivityList,T.Points(:,1),T.Points(:,2),y1)

figure(2)
trisurf(T.ConnectivityList,T.Points(:,1),T.Points(:,2),y2)
zlim ([0 0.8])

[Torsten]

I still wonder what your question has to do with a conditional distribution.

Up to now, you only evaluate two functions (which might represent profit) that depend on two different prices.

But nothing has been said on how probable these prices are.

[Yuriy]

Thank you! I will create a new questions for 3d plots. Regarding probability - I model a competion in 2 periods, where two firms simultanesously draw prices from distribution on interval ($\underbar{p}$ , ) in the first period. The probability to draw a price which is lower than rival's price is . Then the cheapest firm have the following options (whichever is better) in the second period:

1. Increase the price up to rival's price and get a profit y1;
2. Increse the price up to and get a profit y2.

All probable values of these two options I plotted here with your help. And it's only a part of my analysis.

The firm which draw higher price in the first period, always increase the price up to in the second period, however if the cheapest firm chooses option y2, the expensive firm gets additional profits.

So, I am looking for where firms are indifferent in the first period about its price, this means that such distribution of p will represent NE. For that I need to understand how will look and what the lower and upper bound it will have. So far, I found that there are some point on the interval which I call , which affects firms pricing strategies. But I am still looking for NE price distribution which will represent mixed strategies (if it exists).

[Torsten]

The 3d plot is complete (see code above).

[Yuriy]

Thank you!

Is it possible to plot y1 and y2 on the same graph to see how this surfaces interact (cross each other)? Same as it was on 2D plot before where I could see which one is dominating.