Skip to content

Instantly share code, notes, and snippets.

@mbstacy

mbstacy/SOM_MCMC3.m

Created September 10, 2013 20:28
Show Gist options
  • Save mbstacy/6515149 to your computer and use it in GitHub Desktop.
Save mbstacy/6515149 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
SOM_MCMC3.m
% This program is modified by XIA XU from the programm written by Tao Xu to study the inverse problem of
% a one/two/three C pool model using MCMC with initial C pool size
% 06/25/2013-Windows version
clear all;
close all;
format long e;
% random seed is the clock, clock is used because the time we start something is relatively random
RandSeed=clock;
rand('seed',RandSeed(6)); % seed is the first number to start the generator
%Nt = 36; % number of data points
HS = 1000; % histogram
nsim = 5000; % number of simulations
%ctotal = 37.1; % initial C pool size, in mg C g-1 dry soil
% read in data from excel sheet
dataGroup=328;
%read in respiration data
[status_resp,sheet_resp]=xlsfinfo('Days&cumulative respiration-328.xls');
[sheet_r,sheet_c]=size(sheet_resp);
Respdata_sep=cell(sheet_r,sheet_c);
for i=1:sheet_c
Respdata_sep{i}=xlsread('Days&cumulative respiration-328.xls',sheet_resp{i});
end
Respdata=cell(1,dataGroup);
count=1;
for i=1:sheet_c
[resp_r,resp_c]=size(Respdata_sep{i});
col_count=1;
if(i==1) % in MAC version, these 4 lines should be activated and
%keep the the first note colume in Excel (data input-respiration)
resp_c=resp_c-1;
col_count=2;
end
for j=1:(resp_c/2)
Respdata{count}=Respdata_sep{i}(2:end,col_count:(col_count+1));
count=count+1;
col_count=col_count+2;
end
end
% read in SOC data
[status_soc,sheet_soc]=xlsfinfo('SOC-328.xls');
[sheet_r,sheet_c]=size(sheet_soc);
Socdata_sep=cell(sheet_r,sheet_c);
for i=1:sheet_c
Socdata_sep{i}=xlsread('SOC-328.xls',sheet_soc{i});
end
Socdata=zeros(1,dataGroup);
count=1;
for i=1:sheet_c
[m,n]=size(Socdata_sep{i});
for j=1:n
Socdata(count)=Socdata_sep{i}(2,j);
count=count+1;
end
end
%****** prior values (c) and parameter ranges (cmin, cmax) ******************
% f1 f2 parameter names fractions of different pools
a = [0.01 0.15]';
amin = [0 0];
amax = [0.1 0.8];
% k1 k2 k3 respiration rates of different pools
b = [0.001 0.0001 0.000001]';
bmin = [0.0001 0.000 0.000000];
bmax = [0.030 0.0010 0.000004];
% f2,1 f3,1 f1,2 f3,2 f1,3 carbon transfer among pools
c = [0.1 0.002 0.2 0.01 0.2]';
cmin = [0.0 0.000 0.0 0.00 0.0];
cmax = [0.7 0.050 0.4 0.02 0.5];
for group=1:dataGroup
clearvars -except HS RandSeed Respdata Socdata a amax amin b bmax bmin c cmax cmin dataGroup nsim group;
ctotal=Socdata(group);
%*********** read in data ***************
[m,n]=size(Respdata{group});
for i=1:m
if(isnan(Respdata{group}(i,1))==1 || Respdata{group}(i,1)==0)
Nt=i-1;
break;
end
Nt=i;
end
data = Respdata{group}(1:Nt,:);
day = data(:,1); % assign columns to values
%temp = data(:,2)+273.15; % in Kelvin
resp = data(:,2); % daily values of cumulative respiration rates, in mg CO2-C g-1 soil
stdev = std(resp)/sqrt(Nt); % SE of the resp
%simulation starts
adif = (amax-amin)';
bdif = (bmax-bmin)';
cdif = (cmax-cmin)';
a_op=amin'+rand*adif;
a_new=zeros(2,1); % create a new array with zeros
b_op=bmin'+rand*bdif;
J_last = 300000 ; % starting value for last
b_new=zeros(3,1); % create a new array with zeros
record_index=1;
c_op=cmin'+rand*cdif;
c_new=zeros(5,1);
DJ=2*stdev.^2;
%Simulation starts
upgraded=0;
R=8.314; % universal gas constant
for simu=1:nsim
counter=simu
upgraded % shows the upgraded number
%%%%%%%c_new=Generate_2_k(c_op,cmin,cmax); %generate a new point
a_new=Generate_a (a_op,amin, amax);
b_new=Generate_b (b_op,bmin, bmax);
c_new=Generate_c (c_op,cmin, cmax);
f1 = a_new(1);
f2 = a_new(2);
f3 = 1-f1-f2; % not a parameter to estimate: ratio of the least active pool: total pool - ratios of the more active pools
k1 = b_new(1);
k2 = b_new(2);
k3 = b_new(3);
f21 = c_new(1); % carbon transfered from pool 1 to pool 2
f31 = c_new(2); % carbon transfered from pool 1 to pool 3
f12 = c_new(3); % carbon transfered from pool 2 to pool 1
f32 = c_new(4); % carbon transfered from pool 2 to pool 3
f13 = c_new(5); % carbon transfered from pool 3 to pool 1
% model
resp_mod = zeros(day(Nt),1);
CC_mod = zeros(Nt,1);
for i = 1
P1(i) = f1*ctotal-k1*f1*ctotal+f12*k2*f2*ctotal+f13*k3*f3*ctotal;
P2(i) = f2*ctotal-k2*f2*ctotal+f21*k1*f1*ctotal;
P3(i) = f3*ctotal-k3*f3*ctotal+f31*k1*f1*ctotal+f32*k2*f2*ctotal;
R1(i) = (1-f21-f31)*k1*f1*ctotal; %respiration rate of pool1
R2(i) = (1-f12-f32)*k2*f2*ctotal; %respiration rate of pool2
R3(i) = (1-f13)*k3*f3*ctotal; %respiration rate of pool3
C1(i) = R1(i);
C2(i) = R2(i);
C3(i) = R3(i);
fR1(i) = R1(i)/(R1(i)+R2(i)+R3(i)); % proportion of R1 in Rtotal
fR2(i) = R2(i)/(R1(i)+R2(i)+R3(i)); % proportion of R2 in Rtotal
fR3(i) = R3(i)/(R1(i)+R2(i)+R3(i)); % proportion of R3 in Rtotal
resp_mod(i) = C1(i)+C2(i)+C3(i);
end
for i = 2:day(Nt)
P1(i) = P1(i-1)-k1*P1(i-1)+f12*k2*P2(i-1)+f13*k3*P3(i-1); % size of pool1
P2(i) = P2(i-1)-k2*P2(i-1)+f21*k1*P1(i-1); % size of pool2
P3(i) = P3(i-1)-k3*P3(i-1)+f31*k1*P1(i-1)+f32*k2*P2(i-1); % size of pool3
R1(i) = (1-f21-f31)*k1*P1(i-1); %respiration rate of pool1
R2(i) = (1-f12-f32)*k2*P2(i-1); %respiration rate of pool2
R3(i) = (1-f13)*k3*P3(i-1); %respiration rate of pool3
C1(i) = C1(i-1)+R1(i); %cummulative carbon respired from pool1
C2(i) = C2(i-1)+R2(i); %cummulative carbon respired from pool2
C3(i) = C3(i-1)+R3(i); %cummulative carbon respired from pool3
fR1(i) = R1(i)/(R1(i)+R2(i)+R3(i)); % proportion of R1 in Rtotal
fR2(i) = R2(i)/(R1(i)+R2(i)+R3(i)); % proportion of R2 in Rtotal
fR3(i) = R3(i)/(R1(i)+R2(i)+R3(i)); % proportion of R3 in Rtotal
resp_mod(i) = C1(i)+C2(i)+C3(i);
end
for m = 1:Nt
i = day(m);
p1(m)=P1(i);
p2(m)=P2(i);
p3(m)=P3(i);
r1(m)=R1(i);
r2(m)=R2(i);
r3(m)=R3(i);
c1(m)=C1(i);
c2(m)=C2(i);
c3(m)=C3(i);
fr1(m)=fR1(i);
fr2(m)=fR2(i);
fr3(m)=fR3(i);
CC_mod(m) = resp_mod(i);
end
% calculate cost function
J = (norm(CC_mod - resp))^2; % data model error
J_new= J/DJ;
% M-H algorithm
delta_J = J_new-J_last;
%if min(1, exp(-delta_J)) >rand
a_op=a_new;
b_op=b_new;
c_op=c_new;
J_last=J_new;
upgraded=upgraded+1;
a_upgraded(:,upgraded)=a_op;
b_upgraded(:,upgraded)=b_op;
c_upgraded(:,upgraded)=c_op;
CC_record(:,upgraded)=CC_mod;
P1_record(:,upgraded)=p1;
P2_record(:,upgraded)=p2;
P3_record(:,upgraded)=p3;
R1_record(:,upgraded)=r1;
R2_record(:,upgraded)=r2;
R3_record(:,upgraded)=r3;
C1_record(:,upgraded)=c1;
C2_record(:,upgraded)=c2;
C3_record(:,upgraded)=c3;
fR1_record(:,upgraded)=fr1;
fR2_record(:,upgraded)=fr2;
fR3_record(:,upgraded)=fr3;
J_record(:,upgraded) = J_last;
%end
end
%determine if outcomes are accepatable
% parameter ranges
f1=a_upgraded(1,:);
f2=a_upgraded(2,:);
k1=b_upgraded(1,:);
k2=b_upgraded(2,:);
k3=b_upgraded(3,:);
f21=c_upgraded(1,:);
f31=c_upgraded(2,:);
f12=c_upgraded(3,:);
f32=c_upgraded(4,:);
f13=c_upgraded(5,:);
len=length(f1);
a_1=f1(5:len);
a_2=f2(5:len);
b_1=k1(5:len);
b_2=k2(5:len);
b_3=k3(5:len);
c_1=f21(5:len);
c_2=f31(5:len);
c_3=f12(5:len);
c_4=f32(5:len);
c_5=f13(5:len);
groups=50
meanf1 = mean(a_1)
medianf1 = median(a_1)
if (meanf1 - medianf1)< 1*(max(a_1)-min(a_1))/groups
%%%%%%%%%%%%%%%
%%% figures %%%
%%%%%%%%%%%%%%%
F1=figure(1); plot(J_record); %
saveas(F1,strcat('Fig 1_',int2str(group)),'tif');
close(F1);
F2=figure(2);
% plot r1+r2+r3, r_m and r_o
C1_m = mean(C1_record'); % mean modeled respiration rate for the labile pool
C1_m = reshape(C1_m,1,Nt);
C2_m = mean(C2_record');
C2_m = reshape(C2_m,1,Nt);
C3_m = mean(C3_record');
C3_m = reshape(C3_m,1,Nt);
Ctot_m = sum([C1_m(1,:);C2_m(1,:);C3_m(1,:)]); % total modeled respiration, equals resp_m
day_d = reshape(day,1,Nt);
resp_o = reshape(resp,1,Nt);
%temp_d = reshape(temp,1,Nt);
plot(day_d(1,:),C1_m(1,:),'b');hold on;
plot(day_d(1,:),C2_m(1,:),'r');hold on;
plot(day_d(1,:),C3_m(1,:),'g');hold on;
plot(day_d(1,:),Ctot_m(1,:),'c');hold on;
plot(day_d(1,:),resp_o(1,:),'k')
xlabel('time');ylabel('Cummulative C respired (mg C g^-^1 dry soil)'); title(['000' char(176) 'C']);% char is for the degree symbol
legend('C1','C2','C3','Cummu C_m','Cummu C_o');
saveas(F2,strcat('Fig 2_',int2str(group)),'tif');
close(F2);
F3=figure(3);
% parameter ranges
%%f1=a_upgraded(1,:);
%f2=a_upgraded(2,:);
%%k1=b_upgraded(1,:);
%k2=b_upgraded(2,:);
%k3=b_upgraded(3,:);
%f21=c_upgraded(1,:);
%f31=c_upgraded(2,:);
%f12=c_upgraded(3,:);
%f32=c_upgraded(4,:);
%f13=c_upgraded(5,:);
%len=length(f1);
%a_1=f1(5:len);
%a_2=f2(5:len);
%b_1=k1(5:len);
%b_2=k2(5:len);
%b_3=k3(5:len);
%c_1=f21(5:len);
%c_2=f31(5:len);
%c_3=f12(5:len);
%c_4=f32(5:len);
%c_5=f13(5:len);
figure(3);subplot(4,3,1);hist(a_1,30);title('f1');xlabel('Parameter range');ylabel('Frequency');axis([amin(1) amax(1) 1 HS]);
figure(3);subplot(4,3,2);hist(a_2,30);title('f2');xlabel('Parameter range');ylabel('Frequency');axis([amin(2) amax(2) 1 HS]);
figure(3);subplot(4,3,3);hist(b_1,30);title('k1');xlabel('Parameter range');ylabel('Frequency');axis([bmin(1) bmax(1) 1 HS]);
figure(3);subplot(4,3,4);hist(b_2,30);title('k2');xlabel('Parameter range');ylabel('Frequency');axis([bmin(2) bmax(2) 1 HS]);
figure(3);subplot(4,3,5);hist(b_3,30);title('k3');xlabel('Parameter range');ylabel('Frequency');axis([bmin(3) bmax(3) 1 HS]);
figure(3);subplot(4,3,7);hist(c_1,30);title('f2,1');xlabel('Parameter range');ylabel('Frequency');axis([cmin(1) cmax(1) 1 HS]);
figure(3);subplot(4,3,8);hist(c_2,30);title('f3,1');xlabel('Parameter range');ylabel('Frequency');axis([cmin(2) cmax(2) 1 HS]);
figure(3);subplot(4,3,9);hist(c_3,30);title('f1,2');xlabel('Parameter range');ylabel('Frequency');axis([cmin(3) cmax(3) 1 HS]);
figure(3);subplot(4,3,10);hist(c_4,30);title('f3,2');xlabel('Parameter range');ylabel('Frequency');axis([cmin(4) cmax(4) 1 HS]);
figure(3);subplot(4,3,11);hist(c_5,30);title('f1,3');xlabel('Parameter range');ylabel('Frequency');axis([cmin(5) cmax(5) 1 HS]);
saveas(F3,strcat('Fig 3_',int2str(group)),'tif');
close(F3);
%Mark Stacy Added 9/4/2013
%%groups=50
%meanf1 = mean(a_1)
%medianf1 = median(a_1)
%if (meanf1 - medianf1)< 5*(max(a_1)-min(a_1))/groups
% disp('true')
% disp(meanf1)
% disp(medianf1)
% pause
%end
%pause
%************************
F4=figure(4);
%%% scatter plot of observed soil respiration vs modeled soil respiration
C_m = mean(CC_record');
C_m = reshape(C_m,1,Nt);
plot(resp_o(1,:),C_m(1,:),'.');xlabel('Observed cummulative C');ylabel('Modeled cummulative C');
refline; %axis([0 20 0 20]);
cor_om = corrcoef(resp_o,C_m) % correlation coeffecient of observed vs modeled resp
x=resp_o(1,:); y=C_m(1,:); % define x and y
R=corrcoef(x,y); R_squared=R(2)^2; % r square in the figure
text(0.2*max(resp_o(1,:)), 0.8*max(C_m(1,:)), ['R^2 = ' num2str(R_squared)]);
saveas(F4,strcat('Fig 4_',int2str(group)),'tif');
close(F4);
F5=figure(5);
% Carbon pool dynamics
P1_m = mean(P1_record'); % size of pool
P1_m = reshape(P1_m,1,Nt); %
P2_m = mean(P2_record'); % size of pool2
P2_m = reshape(P2_m,1,Nt); %
P3_m = mean(P3_record'); % size of pool 3
P3_m = reshape(P3_m,1,Nt); %
Ptotal_m = sum([P1_m(1,:);P2_m(1,:);P3_m(1,:)]); % total C-pool dynamics
plot(day_d(1,:), P1_m(1,:),'b');hold on; xlabel('days');ylabel('C pool dynamics mg C g-1 dry soil');
plot(day_d(1,:), P2_m(1,:),'r'); hold on;
plot(day_d(1,:), P3_m(1,:),'g'); hold on;
plot(day_d(1,:), Ptotal_m(1,:),'c');
legend('P1','P2','P3','Ptot');
saveas(F5,strcat('Fig 5_',int2str(group)),'tif');
close(F5);
F6=figure(6);
R1_m = mean(R1_record'); % respiration rate of pool 1
R1_m = reshape(R1_m,1,Nt); %
R2_m = mean(R2_record'); % respiration rate of pool 2
R2_m = reshape(R2_m,1,Nt); %
R3_m = mean(R3_record'); % respiration rate of pool 3
R3_m = reshape(R3_m,1,Nt); %
Rtotal_m = sum([R1_m(1,:);R2_m(1,:);R3_m(1,:)]);
fR1_m = mean(fR1_record'); % proportional respiration of pool 1
fR1_m = reshape(fR1_m,1,Nt); %
fR2_m = mean(fR2_record'); % proportional respiration of pool 2
fR2_m = reshape(fR2_m,1,Nt); %
fR3_m = mean(fR3_record'); % proportional respiration of pool 3
fR3_m = reshape(fR3_m,1,Nt); %
fRtotal_m = sum([fR1_m(1,:);fR2_m(1,:);fR3_m(1,:)]);
figure(6);subplot(1,2,1);plot(day_d(1,:),R1_m(1,:),'b');hold on;xlabel('days');ylabel('Respiration rate');
figure(6);subplot(1,2,1);plot(day_d(1,:),R2_m(1,:),'r');hold on;
figure(6);subplot(1,2,1);plot(day_d(1,:),R3_m(1,:),'g');hold on;
figure(6);subplot(1,2,1);plot(day_d(1,:),Rtotal_m(1,:),'c');hold on;
legend('R1','R2','R3','Rtot');
figure(6);subplot(1,2,2);plot(day_d(1,:),fR1_m(1,:),'b');hold on;
figure(6);subplot(1,2,2);plot(day_d(1,:),fR2_m(1,:),'r');hold on;
figure(6);subplot(1,2,2);plot(day_d(1,:),fR3_m(1,:),'g');hold on;
figure(6);subplot(1,2,2);plot(day_d(1,:),fRtotal_m(1,:),'c');hold on;
legend('fR1','fR2','fR3','fRtot');
saveas(F6,strcat('Fig 6_',int2str(group)),'tif');
close(F6);
%%%%%%%%%%%%%%%%%%%%%%
% correlations %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%
% correlation coefficients:
% cor12 = corrcoef(c_1,c_2);
% cor13 = corrcoef(c_1,c_3);
% cor23 = corrcoef(c_2,c_3);
%figure(6);
% plot correlations:
% subplot(1,3,1);plot((c_1),(c_2),'.');xlabel('f1');ylabel('k1');title('parameter correlations');
% subplot(1,3,2);plot((c_1),(c_3),'.');xlabel('f1');ylabel('k2');
% subplot(1,3,3);plot((c_2),(c_3),'.');xlabel('k1');ylabel('k2');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% write the parameters and new data into csv files %%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % write 3p_param parameters into csv files
parameters(:,1)=a_1'; % assigns the parameters to one Matrix
parameters(:,2)=a_2';
parameters(:,3)=b_1';
parameters(:,4)=b_2';
parameters(:,5)=b_3';
parameters(:,6)=c_1';
parameters(:,7)=c_2';
parameters(:,8)=c_3';
parameters(:,9)=c_4';
parameters(:,10)=c_5';
colnames = {'f1','f2','k1','k2','k3','f21','f31','f12','f32','f13'};
fid1=fopen(strcat('3p_param_',int2str(group),'.csv'),'wt');
fprintf(fid1,'%s,%s,%s,%s,%s,%s,%s,%s,%s,%s\n','f1','f2','k1','k2','k3','f21','f31','f12','f32','f13');
[rows,cols]=size(parameters);
for i=1:rows
for j=1:cols
if(j<cols)
fprintf(fid1,'%f,',parameters(i,j));
else
fprintf(fid1,'%f\n',parameters(i,j));
end
end
%fprintf(fid1,'%f\n',parameters(i,j));
end
fclose(fid1);
% write resp_mod to csv
dat_m(:,1)=day_d(1,:);
%dat_m(:,2)= temp_d(1,:)-273.15;
dat_m(:,2)=P1_m(1,:);
dat_m(:,3)=P2_m(1,:);
dat_m(:,4)=P3_m(1,:);
dat_m(:,5)=R1_m(1,:);
dat_m(:,6)=R2_m(1,:);
dat_m(:,7)=R3_m(1,:);
dat_m(:,8)=fR1_m(1,:);
dat_m(:,9)=fR2_m(1,:);
dat_m(:,10)=fR3_m(1,:);
dat_m(:,11)=C1_m(1,:);
dat_m(:,12)=C2_m(1,:);
dat_m(:,13)=C3_m(1,:);
dat_m(:,14)=Ctot_m(1,:);
dat_m(:,15)=resp_o(1,:);
%colnames2 = {'day','temp','P1_m','P2_m','P3_m','R1_m','R2_m','R3_m','fR1_m','fR2_m','fR3_m','C1_m','C2_m','C3_m','Ctot_m','Ctot_o'};
%xlswrite('3p_resp_m.csv',colnames2,1,'A1');
%xlswrite('3p_resp_m.csv',dat_m,1,'A2')
fid2=fopen(strcat('3p_resp_m_',int2str(group),'.csv'),'wt');
fprintf(fid2,'%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s\n','day','P1_m','P2_m','P3_m','R1_m','R2_m','R3_m','fR1_m',...
'fR2_m','fR3_m','C1_m','C2_m','C3_m','Ctot_m','Ctot_o');
[rows,cols]=size(dat_m);
for i=1:rows
for j=1:cols
if(j<cols)
fprintf(fid2,'%f,',dat_m(i,j));
else
fprintf(fid2,'%f\n',dat_m(i,j));
end
end
%fprintf(fid2,'%f\n',data_m(i,j));
end
fclose(fid2);
% write mle and avg to csv
[fre_a1,val_a1]=hist(a_1,30);
[fre_a2,val_a2]=hist(a_2,30);
[fre_b1,val_b1]=hist(b_1,30);
[fre_b2,val_b2]=hist(b_2,30);
[fre_b3,val_b3]=hist(b_3,30);
[fre_c1,val_c1]=hist(c_1,30);
[fre_c2,val_c2]=hist(c_2,30);
[fre_c3,val_c3]=hist(c_3,30);
[fre_c4,val_c4]=hist(c_4,30);
[fre_c5,val_c5]=hist(c_5,30);
max_fre=0;
max_row=0;
for i=1:length(fre_a1)
if(max_fre<fre_a1(i))
max_fre=fre_a1(i);
max_row=i;
end
end
mu_a1=val_a1(max_row);
max_fre=0;
max_row=0;
for i=1:length(fre_a2)
if(max_fre<fre_a2(i))
max_fre=fre_a2(i);
max_row=i;
end
end
mu_a2=val_a2(max_row);
max_fre=0;
max_row=0;
for i=1:length(fre_b1)
if(max_fre<fre_b1(i))
max_fre=fre_b1(i);
max_row=i;
end
end
mu_b1=val_b1(max_row);
max_fre=0;
max_row=0;
for i=1:length(fre_b2)
if(max_fre<fre_b2(i))
max_fre=fre_b2(i);
max_row=i;
end
end
mu_b2=val_b2(max_row);
max_fre=0;
max_row=0;
for i=1:length(fre_b3)
if(max_fre<fre_b3(i))
max_fre=fre_b3(i);
max_row=i;
end
end
mu_b3=val_b3(max_row);
max_fre=0;
max_row=0;
for i=1:length(fre_c1)
if(max_fre<fre_c1(i))
max_fre=fre_c1(i);
max_row=i;
end
end
mu_c1=val_c1(max_row);
max_fre=0;
max_row=0;
for i=1:length(fre_c2)
if(max_fre<fre_c2(i))
max_fre=fre_c2(i);
max_row=i;
end
end
mu_c2=val_c2(max_row);
max_fre=0;
max_row=0;
for i=1:length(fre_c3)
if(max_fre<fre_c3(i))
max_fre=fre_c3(i);
max_row=i;
end
end
mu_c3=val_c3(max_row);
max_fre=0;
max_row=0;
for i=1:length(fre_c4)
if(max_fre<fre_c4(i))
max_fre=fre_c4(i);
max_row=i;
end
end
mu_c4=val_c4(max_row);
max_fre=0;
max_row=0;
for i=1:length(fre_c5)
if(max_fre<fre_c5(i))
max_fre=fre_c5(i);
max_row=i;
end
end
mu_c5=val_c5(max_row);
avg_a1=mean(a_1);
avg_a2=mean(a_2);
avg_b1=mean(b_1);
avg_b2=mean(b_2);
avg_b3=mean(b_3);
avg_c1=mean(c_1);
avg_c2=mean(c_2);
avg_c3=mean(c_3);
avg_c4=mean(c_4);
avg_c5=mean(c_5);
fid3=fopen(strcat('3p_mle_avg_',int2str(group),'.csv'),'wt');
fprintf(fid3,'%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s\n','mle_f1','mle_f2','mle_k1','mle_k2','mle_k3','mle_f21','mle_f31','mle_f12','mle_f32','mle_f13',...
'avg_f1','avg_f2','avg_k1','avg_k2','avg_k3','avg_f21','avg_f31','avg_f12','avg_f32','avg_f13');
fprintf(fid3,'%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f\n',mu_a1,mu_a2,mu_b1,mu_b2,mu_b3,mu_c1,mu_c2,mu_c3,mu_c4,mu_c5,...
avg_a1,avg_a2,avg_b1,avg_b2,avg_b3,avg_c1,avg_c2,avg_c3,avg_c4,avg_c5);
fclose(fid3);
end %end of output if
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment