How to loop function 1000 times

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

Answer 1

[Kavya Vuriti]

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.