%PROGRAM NAME: biochar.m
%Biochar data set with missing values =0): 
% data(j,k) data(j,k+1) represent flux (slope) for jth week, for kth time,
% and for replications 1 and 2.
clear all
data=[1075.734858	1785.020358	    2079.875029	  2577.519237	2860.690076	  3452.196667	3443.486879	 4095.481276;
      3619.414529	3881.481338		4164.030485	  4596.624326	4856.147019	  5058.525969	5610.318515  0
      2746.440677	2545.395841	    3455.029268	  3324.402792	4126.142597	  3953.044998	5052.825716	 4652.149714;
      3399.420053	3275.072076	    3700.210581	  3752.381658	4441.751425	  4330.288838	4549.429985	 4828.779046;
	  3459.368753	3842.1574		4778.570178	  4561.973148   0             0             0            0;
      3823.689766	3622.763143	    4518.898143	  4603.357772	4896.995311	  5190.942931	5426.633452	 5344.129114;
      5292.25713	4475.161472	    5352.887028	  4980.970289	5697.803657	  5194.280558	5817.607871	 5362.814048;
      6189.584429	6255.531255	    6001.092726	  7349.396047	6902.202019	  8281.503115	7658.594124  0];
  %Modify data array to fill in missing values by the averave of the
  %corresponding column:
  data(2,8)=sum(data(:,8))/5; data(5,8)=data(2,8); data(8,8)=data(2,8);
  for k=5:7
      data(5,k)=sum(data(:,k))/7;
  end
  Data=data'; t=[0 0 15 15 30 30 45 45]'; tR=t(1:2:8);
  a=zeros(1,8);   b=zeros(1,8);
  Chat=zeros(2,2,8); Cabhat=Chat;
  for k=1:8  %Compute linear models for each week
      D=[t,Data(:,k)];
      Chat(:,:,k)=cov(D,1);
      a(k)=Chat(1,2,k)/Chat(1,1,k); %Estimated slope for kth week
      M=mean(D);
      b(k)=M(2)-a(k)*M(1); %Estimated intercept for kth week
      %-------------------------------------------------------
      %Compute [ahat,bhat] covariance matrix:
      Cabhat(1,1,k)=Chat(2,2,k)/(8*Chat(1,1,k));
      Cabhat(2,2,k)=Cabhat(1,1,k)*mean(t.^2);
      Cabhat(1,2,k)=Cabhat(1,1,k)*(-mean(t));
      Cabhat(2,1,k)=Cabhat(1,2,k);
  end
%====================================
a=a'; b=b';
 Bhat =kron(a,tR')+kron(b,ones(1,4)); 
 figure(1)
 plot(tR',Bhat)
 grid
 title('Biochar Prediction Models for Each of the 8 Weeks')
 xlabel('Time (minutes)')
 M=num2str(a);
 legend(M,'Location','SouthEast')
 ylabel('Biochar (ppm)')
 pause
 figure(2)
 for k=1:8
     k
     vBhat=tR'.^2*Cabhat(1,1,k)+Cabhat(2,2,k)+2*tR'*Cabhat(1,2,k);
     std2Bhat=2*sqrt(vBhat);
     BB=Bhat(k,:);
     B1=BB-std2Bhat;
     B2=BB+std2Bhat;
     if k>1
         hold on
     end
     plot(tR',[BB;B1;B2])
     title('Biochar +/- 2sigma Prediction Bounds for kth week')
     xlabel('Time (minuts)')
     ylabel('Biochar (ppm)')
     pause
 end
 %Investigate Prediction Errors:
 BBhat=[Bhat,Bhat];
 Q=data - BBhat;
 figure(3)
 for m=1:2
     m1=1+(m-1)*4;
     m2=m1+3;
 subplot(2,1,m),plot(tR',Q(1:2:8,m1:m2),'*-')
 hold on
  plot(tR',Q(2:2:8,1:4),'*--')
 grid
 title(['Biochar Replicate Model Errors for Each Weeks ',num2str(m1),' to ', num2str(m2),' '])
 xlabel('Time (minutes)')
 W=m1:m2; W=[W W];
 W=num2str(W');
 legend(W,'Location','EastOutside')
 ylabel('Biochar Error (ppm)')
 end
 pause
 %-----------------------
 Ma = mean(a)
 STDa=std(a)
 figure(4)
 hist(a)
 title('Histogram of ahat for n=8 weeks')
 grid
 pause
 figure(5)
 plot(a,'*-')
 title('Plot of ahat(n) for n=1:8 weeks')
 xlabel('week')
 grid
 
 
  
  
