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
