CODEWITHSUNDEEP
matlabmatlab-figure
David
David

Reputation: 77

compute the average over simulation for different parameters values

I am working on something similar to the following example: I want to compute <x(t)> that is the average of a function x(t) over number of simulations. To do this, I generate the following code:

sim=50;% number of simulations 
t=linspace(0,1);% time interval

a_range=[1,2,3];% different values for the parameter a
b_range=[0,0.5,1];% different values for the parameter b
z=zeros(1,sim);
theta=zeros(1,sim);

for nplot=1:3
   a=a_range(nplot);
   b=b_range(nplot);
   average_x=zeros(nplot,sim);
   for i=1:sim
       z(i)=rand(1);% random number for every simulation
       theta(i)=pi*rand(1);% random number for every simulation    
    x=z(i)*t.^2+a*sin(theta(i))+b.*tan(theta(i));% the function

   end
   average_x(nplot,sim)=mean(x);% average over the number of simulations
end

fname=['xsin.mat'];
save(fname)

The time is a vector 1 by 100 and x is a vector 1 by 100, and average_x is 1 by 50. What I am looking for is to write a script to load the file and plot the average against time for different parameters a and b. So I want to write a code to generate three figures such that in figure 1 I will plot the average

plot(t,average_x)

for a=1 and b=0.

Then in figure 2 I will plot the average again but for a=2 and b=0.5 and so on. The problem is the dimension of time t and the average are not the same. How can I fix this problem and generate three distinct figures.

Upvotes: 0

Views: 218

Answers (1)

EBH
EBH

Reputation: 10440

If I got your intention correctly, this is what you look for:

sim = 50;% number of simulations 
t = linspace(0,1);% time interval

a_range = [1,2,3];% different values for the parameter a
b_range = [0,0.5,1];% different values for the parameter b

% NO NEED TO GENERATE THE RANDOM NUMBERS ONE BY ONE:
theta = pi*rand(sim,1);% random number for every simulation
z = rand(sim,1); % random number for every simulation

% YOU SOULD INITIALIZE ALL YOUR VARIABLES OUTSIDE THE LOOPS:
x = zeros(sim,numel(t)); 
average_x = zeros(3,numel(t));% the mean accross simulations
% for average accros time use:
% average_x = zeros(3,sim);

for nplot=1:3
   a = a_range(nplot);
   b = b_range(nplot);
   for i=1:sim
       x(i,:) = z(i)*t.^2+a*sin(theta(i))+b.*tan(theta(i));% the function
   end
   average_x(nplot,:) = mean(x); % average over the number of simulations
   % average_x(nplot,:) = mean(x,2); % average accross time
end
% save the relevant variables:
save('results.mat','average_x','t')

In another file you can write:

load('results.mat')
for k = 1:size(average_x,1)
    figure(k)
    plot(t,average_x(k,:))
    title(['Parameter set ' num2str(k)])
    xlabel('Time')
    ylabel('mean x')
end

This is the plot in one figure (if you want then average over simulations):

mean sim

BTW, if you want to make your code more compact and fast, you can vectorize it, mainly using bsxfun. Here is a demonstration with your code:

% assuming all parameters are defined as above:
zt = bsxfun(@times,z,t.^2); % first part of the function 'z(i)*t.^2'
% second part of the function 'a*sin(theta(i)) + b.*tan(theta(i))':
ab = bsxfun(@times,a_range,sin(theta)) + bsxfun(@times,b_range,tan(theta));
% convert the second part to the right dimensions and size:
ab = repmat(reshape(ab,[],1,3),1,numel(t),1); 
x = bsxfun(@plus,zt,ab); % the function
average_x = squeeze(mean(x)); % take the mean by simulation
plot(t,average_x) % plot it all at once, as in the figure above
xlabel('Time')
ylabel('mean x')

Upvotes: 2

Related Questions