3-d plot generated from simulations

Hi all,

Could someone please tell me how to generate a graph like the one in the figure below?

My idea is that I have a variable X that I simulate 1000 times with a normal distribution to which I want to vary the mean (from 75 to 125) and the standard deviation (from 1 to 50). From the different values of X obtained, I would like to calculate the Y variable. My idea is to make a graph where the x-axis is the different mean values, the y-axis is the different standard deviation values from X, and the Z-axis is the expected value (Z) obtained for the Y function.

n=1000; %number of sample
X = normrnd(100,40,n,1); %Demand Simulation,here I want to change the mean and the std. deviation
Y=0.3*((X/10).^(2)); 
Z=mean(Y);

Thank you very much for your help!

[Scott MacKenzie]

I see a couple of issues in your question. First, it seems Z is only a functon of X. Second, X, in your code, is 1000 values from a normal distribution, but, X, in your question varies from 75 to 125. Can you clarify?

[Angelavtc]

Of course @Scott MacKenzie. Indeed, Z is only a function of X. And in my example I simulate 1000 values of X but keep the mean and standard deviation fixed (at 100 and 40 respectively). What I need is for these two values to vary, the first between 75 and 125 and the second between 1 and 50. I imagine this is accomplished with a for, but I can't think of how to link it with a plot. Indeed, the simulation of the 1000 values is incorrect for the purposes of my plot because I forgot to add that, from these 1000 simulations of X going into Z what I am interested in plotting is the expected value of Z. I am correcting this now. Thank you!

[Scott MacKenzie]

I'm still scratching my head over your question. BTW, your code reduces to...

n=1000; %number of sample
X = normrnd(100,40,n,1); %Demand Simulation,here I want to change the mean and the std. deviation
Z=mean(0.3*((X/10).^(2))); 

I don't see the connection between this code (with two variables) and the 3D plot you are trying to create (with three variables).

The plot you are trying to create shows z as function of x (which varies from 75 to 125 -- evidently a mean) and y (which varies from 1 to 50 -- evidently a standard deviation). It seems z is a function of two values, one a mean and the other a standard deviation. If that's the case, why are you using normrnd and generate 1000 sample points from the distribution. It's only the mean and standard deviation that you need, and these are known.

[Angelavtc]

Dear @Scott MacKenzie...Each observation of 1000 points generated by x represents for me a demand value where (in the case of normrnd(100,40,n,1)) together they have a mean of 100 and a standard deviation of 40. Thus, to this simulation of 1000 points correspond 1000 different prices (my Y) and there is a specific mean for these 1000 prices (what I call Z). I am interested to know what happens to that mean in the price (my Z) when I change my demand simulation to another one with 1000 different observations, for example, with a mean of 100 and a standard deviation of 41and so on until all possible combinations of mean (75 to 125) vs standard deviation are exhausted (1 to 50). I don't quite understand what you mean? if you could explain it better, I would appreciate it. Thank you!

[Scott MacKenzie]

Let's try this. In your example plot, one of the x-axis values in 99 and one of the y-axis values is 28. It's hard to tell, but z looks to be about -1.5 or so. Can you give the formula for z in terms of x and y?

[Angelavtc]

@Scott MacKenzie ok I think I'm close after all your hints... it's something like this? and how to graph it?

iMax = 50;
jMax = 50;
myArray = zeros(iMax, jMax);
for iIdx = 1:iMax
    for jIdx = 1:jMax
       myArray(iIdx,jIdx) =mean(0.3*((normrnd(iIdx,jIdx,n,1)/10).^(2))); 
    end
end

[Scott MacKenzie]

OK, we're getting there. I think you forgot n=1000; at the beginning of your code. If you add that and the following lines at the end...

[X, Y] = meshgrid(1:iMax, 1:jMax);
surf(X, Y, myArray');
xlabel('X'); ylabel('Y'); zlabel('Z');

here's the resulting plot:

This gets you closer to your goal. The Z values are no where near those in your example plot and the surface is quite rough. The roughness is perhaps an expected side effect of using normrnd.

[Angelavtc]

Thank you very much @Scott MacKenzie! this works very well. The graph is not the same as the one I put in the example, as I only used it as a reference for my own function. One last question, how do I make mean (iIdx) start at value 75 and end at 125? Because with my example it starts from 1.

[Angelavtc]

Additionally could you please copy and paste your answer so I can accept it? @Scott MacKenzie

[Scott MacKenzie]

OK, sure. I don't this is the exactly what you were looking for, but hopefully it helps. Good luck.

[Angelavtc]

@Scott MacKenzie just could you please explain me how to change the loop in order to make the mean of the normal distribution starts at 75 instead of 1? Cant figure it out :/ Thank you!

[Scott MacKenzie]

Yea, I was fiddling around with that also. Here's what I put together (with the variable names changed to reflect the x, y, and z axes):

n = 1000;
x = 75:125;
y = 1:50;

for i = 1:length(x) % mean
    for j = 1:length(y) % sd
       Z(i,j) = mean(0.3 * ((normrnd(x(i),y(j),n,1) / 10).^(2)));      
    end   
end

[X, Y] = ndgrid(x, y);

h = surf(X, flip(Y), Z);
xlabel('X'); ylabel('Y'); zlabel('Z');

I still haven't figure out how to reverse the order of the y axis values.

[Angelavtc]

Thank you @Scott MacKenzie! this is wondeful :)

[Scott MacKenzie]

You're welcome. But, please be aware that by flipping Y (or using the transpose of myArray in my original answer), the relationship between the data and the changed variable is also flipped. Really, I was just trying to create a reasonable visual graph.

BTW, to smooth the surface, you can add this after the outer for-loop:

Z = smoothdata(Z);

[Angelavtc]

Dear @Scott MacKenzie, thanks for the tip, is this the reason why when I replicate the graph in gray (here below) I get negative results on the z-axis (bias in forward price) instead of positive ones? If so, how can I correct this? I am interested in replicating exactly this graph, with the same order in the axis. I also notice that my demand standard deviation axis is upside down with respect to the original graph, how can I change this? Thank you!

[Scott MacKenzie]

Hmm, not sure. Are you using a different formula? I think there is something fundamentally wrong in the calculations, as opposed an issue in plotting the results. Even the data in my solution (created from the nested for-loop you provided) don't come close to matching the data in the original gray scale graph. So, if you can't recreate the data, there's no change of recreating the plot.

[Angelavtc]

@Scott MacKenzie That's right, I'm using a different set of formulas that I didn't put in the question because they make it complicated and don't add much to my graphical doubt. So maybe it's a problem in my formulas that I'm going to check. However, in just graphical terms do you know how to change the order of the y-axis so that it starts at 40 and ends at 0 (descending order) and could you elaborate more on what you mean by the y-flip problem? Thanks in advance!

[Scott MacKenzie]

@Angelavtc I just dug a little deeper and figure this out. There is an axis property that allows you to reverse the direction of the data along an axis. After the surf function, put...

set(gca, 'YDir', 'reverse');

[Angelavtc]

@Scott MacKenzie the reverse command works very well, thank you. Then flip only inverts the elements... The strange thing is that when I plot with or without flip, the graph does not change.

[Scott MacKenzie]

With using flip on a matrix, you need to be aware of whether you are flipping along the rows or column. flip(Y) or flip(Y,1) flips the rows, but flip(Y,2) flips the columns. Try the latter and you'll see a difference.

[Angelavtc]

Ok! Thanks for everything @Scott MacKenzie I really appreciate all your help and wish you the best :)

Answer 1

[Scott MacKenzie]

n = 1000;
iMax = 50;
jMax = 50;
myArray = zeros(iMax, jMax);
for iIdx = 1:iMax
    for jIdx = 1:jMax
       myArray(iIdx,jIdx) = mean(0.3*((normrnd(iIdx,jIdx,n,1)/10).^(2))); 
    end
end

[X, Y] = ndgrid(1:iMax, 1:jMax);
surf(X, Y, myArray');
xlabel('X'); ylabel('Y'); zlabel('Z');