How to loop function 1000 times
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
Hello,
I have just created this script that determines the weight of phosphorous entering lake Erie each month. It determines the monthly amounts of phosphorous for n = 20 years. I need to loop this script 1000 times in order to get more accurate average values for each month.
Thanks
Here is the script:
clear
clc
mu = [1.606396186,5.090733506,5.738006686,15.84,5.025985475,4.394665359,3.466560461,0.26];
sigma = [0.220905289,1.190553451,0.709631846,6.07,1.125946115,1.346005319,1.533483881,0.12];
bc = [0,1,1,.5,1,1,1,-.5];
A = 21.88;
B = 0.23;
n = 20; % Number of years
for i = 1:n
p = rand(n,8); % Random number generator, p is a n-by-8 matrix
% Data transformation
january_trans(i) = gaminv(p(n,1),A,B);
february_trans(i) = logninv(p(n,2),mu(2),sigma(2));
march_trans(i) = logninv(p(n,3),mu(3),sigma(3));
april_trans(i) = norminv(p(n,4),mu(4),sigma(4));
may_trans(i) = logninv(p(n,5),mu(5),sigma(5));
june_trans(i) = logninv(p(n,6),mu(6),sigma(6));
july_trans(i) = logninv(p(n,7),mu(7),sigma(7));
august_trans(i) = norminv(p(n,8),mu(8),sigma(8));
% Storing transformed data, (row,column)=(month,year)
transformed(:,i) = [january_trans(i),february_trans(i),march_trans(i),april_trans(i),may_trans(i),june_trans(i),july_trans(i),august_trans(i)];
% Data untransformation
january_untrans(i) = exp(january_trans(i));
february_untrans(i) = february_trans(i).^(1/bc(2));
march_untrans(i) = march_trans(i).^(1/bc(3));
april_untrans(i) = april_trans(i).^(1/bc(4));
may_untrans(i) = may_trans(i).^(1/bc(5));
june_untrans(i) = june_trans(i).^(1/bc(6));
july_untrans(i) = july_trans(i).^(1/bc(7));
august_untrans(i) = august_trans(i).^(1/bc(8));
% Storing untransformed data, (row,column)=(month,year)
untransformed(:,i) = [january_untrans(i),february_untrans(i),march_untrans(i),april_untrans(i),may_untrans(i),june_untrans(i),july_untrans(i),august_untrans(i)];
end
yearly_sum = sum(untransformed)';
for j = 1:n
Beta_o = [-39.6,-10.1];
Beta_w = [56.7,113];
Beta_b = [4.21,9.71];
Beta_t = [.77,5.16];
Beta_g = [1.61,3.69];
Sigma_g = [1.63,4.71];
Sigma_e = [.15,5.78];
% Random number generator between [min,max]
beta_o(j) = Beta_o(1)+rand(1)*(Beta_o(2)-Beta_o(1));
beta_w(j) = Beta_w(1)+rand(1)*(Beta_w(2)-Beta_w(1));
beta_b(j) = Beta_b(1)+rand(1)*(Beta_b(2)-Beta_b(1));
beta_t(j) = Beta_t(1)+rand(1)*(Beta_t(2)-Beta_t(1));
beta_g(j) = Beta_g(1)+rand(1)*(Beta_g(2)-Beta_g(1));
sigma_g(j) = Sigma_g(1)+rand(1)*(Sigma_g(2)-Sigma_g(1));
sigma_e(j) = Sigma_e(1)+rand(1)*(Sigma_e(2)-Sigma_e(1));
w_i = untransformed(1:6,:)./1000;
for k = 1:6
m(k) = k;
if (m(k) < beta_g(j)) && (m(k) >= (beta_g(j)-1))
psi(k,j) = m(k) + 1 - beta_g(j);
elseif (m(k) >= beta_g(j))
psi(k,j) = 1;
else
psi(k,j) = 0;
end
end
end
psi_sum = sum(psi);
product = w_i.*psi;
product_sum = sum(product);
for ii = 1:n
W_i(ii) = (1/psi_sum(ii))*product_sum(ii);
end
T_i = 0:1:19;
for jj = 1:n
if (beta_o(jj)+beta_w(jj)*W_i(jj)+beta_t(jj)*T_i(jj) > 0)
z_hat_i(jj) = beta_b(jj)+beta_o(jj)+beta_w(jj)*W_i(jj)+beta_t(jj)*T_i(jj);
else
z_hat_i(jj) = beta_b(jj);
end
A_1(jj) = (z_hat_i(jj)^2)/(sigma_g(jj)^2);
B_1(jj) = (sigma_g(jj)^2)/z_hat_i(jj);
gamma_1(jj) = gaminv(rand(1),A_1(jj),B_1(jj))-z_hat_i(jj);
A_2(jj) = ((z_hat_i(jj)+gamma_1(jj))^2)/(sigma_e(jj)^2);
B_2(jj) = (sigma_e(jj)^2)/(z_hat_i(jj)+gamma_1(jj));
z_i(jj) = gaminv(rand(1),A_2(jj),B_2(jj));
end
0 Kommentare
Antworten (1)
Kavya Vuriti
am 11 Nov. 2019
Hi,
You can try writing a MATLAB function defining the code in your script. Output of the function can be monthly amount of phosphorus for 20 years. Assuming execution time is not a concern, you can write a ‘for’ loop where output of the function is taken and accumulated for each month. These accumulated sums can be used to find average value for each month.
0 Kommentare
Siehe auch
Kategorien
Mehr zu Entering Commands finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!