異なるパラメータ値のシミュレーション上の平均を計算する

Question

私は次の例のようなものに取り組んでいます：<x(t)>は、シミュレーションの数を超えた関数x(t)の平均です。これを行うために、私は次のコードを生成します。

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) 

時間は100で、ベクトル1であり、xは100で、ベクトル1である、とaverage_xは私が探しています何50で1で書くことですスクリプトを使用してファイルをロードし、異なるパラメータaとbの平均時間をプロットします。だから私は図1で平均をプロットするように3つの数字を生成するコードを書いて、a = 1とb = 0に対して平均をプロットしたいと思います。

次に、図2では、平均を再びプロットしますが、a = 2およびb = 0.5などとします。問題は時間tの次元で平均は同じではありません。どのようにしてこの問題を解決し、3つの異なる数字を生成できますか？あなたがしたい場合

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 

これは1つの図のプロットである（：

Answer 1

私が正しくあなたの意図を持っている場合、これはあなたが探して何をされて：あなたが書くことができ、別のファイルでは

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') 

その後、シミュレーションの平均）：

---

ところで、コードをもっとコンパクトにするには、bsxfunを中心にベクトル化することができます。

% 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')