Step 1: Organize the input data in a text file (e.g. data.txt) Create a text file named data.txt copy and paste the two columns of data (one for forecasts, the other for observations) into the data.txt file.
tips: only put in the data values, please don't add any character or string. The data.txt file is uploaded. Step 2: Create MATLAB script and read in dataOpen MATLAB, and follow the steps: File->New->Script Type codes in the new script: %first column, observations; second column, forecast
fp=fopen('data.txt','r');
data=zeros(196,2);
for i = 1 : 196
for j = 1 : 2
data(i,j)=fscanf(fp,'%d',1);%for temperature
end
end
for i = 1 : 196
for j = 1 : 2
data2(i,j)=fscanf(fp,'%d',1);%for precipitation
end
end Step 3: Get contingency tablesTwo contingency tables, for temperature 9x9 and precipitation 2x5 respectively. Add one more row and one more column to store the summation of each column or row. So they are tables of 10x10 and 3x6. Group the data into corresponding slots in the tables. Finally, sum up each row and column. table=zeros(10,10); table2=zeros(3,6); %%contingency table 1 temperature for i = 1 : 196 if (data(i,1)>=10 && data(i,1)<=19) if (data(i,2)>=10 && data(i,2)<=19) table(1,1)=table(1,1)+1; elseif (data(i,2)>=20 && data(i,2)<=29) table(1,2)=table(1,2)+1; elseif (data(i,2)>=30 && data(i,2)<=39) table(1,3)=table(1,3)+1; elseif (data(i,2)>=40 && data(i,2)<=49) table(1,4)=table(1,4)+1; elseif (data(i,2)>=50 && data(i,2)<=59) table(1,5)=table(1,5)+1; elseif (data(i,2)>=60 && data(i,2)<=69) table(1,6)=table(1,6)+1; elseif (data(i,2)>=70 && data(i,2)<=79) table(1,7)=table(1,7)+1; elseif (data(i,2)>=80 && data(i,2)<=89) table(1,8)=table(1,8)+1; elseif (data(i,2)>=90) table(1,9)=table(1,9)+1; end elseif (data(i,1)>=20 && data(i,1)<=29) if (data(i,2)>=10 && data(i,2)<=19) table(2,1)=table(2,1)+1; elseif (data(i,2)>=20 && data(i,2)<=29) table(2,2)=table(2,2)+1; elseif (data(i,2)>=30 && data(i,2)<=39) table(2,3)=table(2,3)+1; elseif (data(i,2)>=40 && data(i,2)<=49) table(2,4)=table(2,4)+1; elseif (data(i,2)>=50 && data(i,2)<=59) table(2,5)=table(2,5)+1; elseif (data(i,2)>=60 && data(i,2)<=69) table(2,6)=table(2,6)+1; elseif (data(i,2)>=70 && data(i,2)<=79) table(2,7)=table(2,7)+1; elseif (data(i,2)>=80 && data(i,2)<=89) table(2,8)=table(2,8)+1; elseif (data(i,2)>=90) table(2,9)=table(2,9)+1; end elseif (data(i,1)>=30 && data(i,1)<=39) if (data(i,2)>=10 && data(i,2)<=19) table(3,1)=table(3,1)+1; elseif (data(i,2)>=20 && data(i,2)<=29) table(3,2)=table(3,2)+1; elseif (data(i,2)>=30 && data(i,2)<=39) table(3,3)=table(3,3)+1; elseif (data(i,2)>=40 && data(i,2)<=49) table(3,4)=table(3,4)+1; elseif (data(i,2)>=50 && data(i,2)<=59) table(3,5)=table(3,5)+1; elseif (data(i,2)>=60 && data(i,2)<=69) table(3,6)=table(3,6)+1; elseif (data(i,2)>=70 && data(i,2)<=79) table(3,7)=table(3,7)+1; elseif (data(i,2)>=80 && data(i,2)<=89) table(3,8)=table(3,8)+1; elseif (data(i,2)>=90) table(3,9)=table(3,9)+1; end elseif (data(i,1)>=40 && data(i,1)<=49) if (data(i,2)>=10 && data(i,2)<=19) table(4,1)=table(4,1)+1; elseif (data(i,2)>=20 && data(i,2)<=29) table(4,2)=table(4,2)+1; elseif (data(i,2)>=30 && data(i,2)<=39) table(4,3)=table(4,3)+1; elseif (data(i,2)>=40 && data(i,2)<=49) table(4,4)=table(4,4)+1; elseif (data(i,2)>=50 && data(i,2)<=59) table(4,5)=table(4,5)+1; elseif (data(i,2)>=60 && data(i,2)<=69) table(4,6)=table(4,6)+1; elseif (data(i,2)>=70 && data(i,2)<=79) table(4,7)=table(4,7)+1; elseif (data(i,2)>=80 && data(i,2)<=89) table(4,8)=table(4,8)+1; elseif (data(i,2)>=90) table(4,9)=table(4,9)+1; end elseif (data(i,1)>=50 && data(i,1)<=59) if (data(i,2)>=10 && data(i,2)<=19) table(5,1)=table(5,1)+1; elseif (data(i,2)>=20 && data(i,2)<=29) table(5,2)=table(5,2)+1; elseif (data(i,2)>=30 && data(i,2)<=39) table(5,3)=table(5,3)+1; elseif (data(i,2)>=40 && data(i,2)<=49) table(5,4)=table(5,4)+1; elseif (data(i,2)>=50 && data(i,2)<=59) table(5,5)=table(5,5)+1; elseif (data(i,2)>=60 && data(i,2)<=69) table(5,6)=table(5,6)+1; elseif (data(i,2)>=70 && data(i,2)<=79) table(5,7)=table(5,7)+1; elseif (data(i,2)>=80 && data(i,2)<=89) table(5,8)=table(5,8)+1; elseif (data(i,2)>=90) table(5,9)=table(5,9)+1; end elseif (data(i,1)>=60 && data(i,1)<=69) if (data(i,2)>=10 && data(i,2)<=19) table(6,1)=table(6,1)+1; elseif (data(i,2)>=20 && data(i,2)<=29) table(6,2)=table(6,2)+1; elseif (data(i,2)>=30 && data(i,2)<=39) table(6,3)=table(6,3)+1; elseif (data(i,2)>=40 && data(i,2)<=49) table(6,4)=table(6,4)+1; elseif (data(i,2)>=50 && data(i,2)<=59) table(6,5)=table(6,5)+1; elseif (data(i,2)>=60 && data(i,2)<=69) table(6,6)=table(6,6)+1; elseif (data(i,2)>=70 && data(i,2)<=79) table(6,7)=table(6,7)+1; elseif (data(i,2)>=80 && data(i,2)<=89) table(6,8)=table(6,8)+1; elseif (data(i,2)>=90) table(6,9)=table(6,9)+1; end elseif (data(i,1)>=70 && data(i,1)<=79) if (data(i,2)>=10 && data(i,2)<=19) table(7,1)=table(7,1)+1; elseif (data(i,2)>=20 && data(i,2)<=29) table(7,2)=table(7,2)+1; elseif (data(i,2)>=30 && data(i,2)<=39) table(7,3)=table(7,3)+1; elseif (data(i,2)>=40 && data(i,2)<=49) table(7,4)=table(7,4)+1; elseif (data(i,2)>=50 && data(i,2)<=59) table(7,5)=table(7,5)+1; elseif (data(i,2)>=60 && data(i,2)<=69) table(7,6)=table(7,6)+1; elseif (data(i,2)>=70 && data(i,2)<=79) table(7,7)=table(7,7)+1; elseif (data(i,2)>=80 && data(i,2)<=89) table(7,8)=table(7,8)+1; elseif (data(i,2)>=90) table(7,9)=table(7,9)+1; end elseif (data(i,1)>=80 && data(i,1)<=89) if (data(i,2)>=10 && data(i,2)<=19) table(8,1)=table(8,1)+1; elseif (data(i,2)>=20 && data(i,2)<=29) table(8,2)=table(8,2)+1; elseif (data(i,2)>=30 && data(i,2)<=39) table(8,3)=table(8,3)+1; elseif (data(i,2)>=40 && data(i,2)<=49) table(8,4)=table(8,4)+1; elseif (data(i,2)>=50 && data(i,2)<=59) table(8,5)=table(8,5)+1; elseif (data(i,2)>=60 && data(i,2)<=69) table(8,6)=table(8,6)+1; elseif (data(i,2)>=70 && data(i,2)<=79) table(8,7)=table(8,7)+1; elseif (data(i,2)>=80 && data(i,2)<=89) table(8,8)=table(8,8)+1; elseif (data(i,2)>=90) table(8,9)=table(8,9)+1; end elseif (data(i,1)>=90) if (data(i,2)>=10 && data(i,2)<=19) table(9,1)=table(9,1)+1; elseif (data(i,2)>=20 && data(i,2)<=29) table(9,2)=table(9,2)+1; elseif (data(i,2)>=30 && data(i,2)<=39) table(9,3)=table(9,3)+1; elseif (data(i,2)>=40 && data(i,2)<=49) table(9,4)=table(9,4)+1; elseif (data(i,2)>=50 && data(i,2)<=59) table(9,5)=table(9,5)+1; elseif (data(i,2)>=60 && data(i,2)<=69) table(9,6)=table(9,6)+1; elseif (data(i,2)>=70 && data(i,2)<=79) table(9,7)=table(9,7)+1; elseif (data(i,2)>=80 && data(i,2)<=89) table(9,8)=table(9,8)+1; elseif (data(i,2)>=90) table(9,9)=table(9,9)+1; end end end %%contingency table 2 precipitation for i = 1 : 196 if (data2(i,1)==0) if (data2(i,2)>=0 && data2(i,2)<=20) table2(1,1)=table2(1,1)+1; elseif (data2(i,2)>=21 && data2(i,2)<=40) table2(1,2)=table2(1,2)+1; elseif (data2(i,2)>=41 && data2(i,2)<=60) table2(1,3)=table2(1,3)+1; elseif (data2(i,2)>=61 && data2(i,2)<=80) table2(1,4)=table2(1,4)+1; elseif (data2(i,2)>=81 && data2(i,2)<=100) table2(1,5)=table2(1,5)+1; end elseif (data2(i,1)==100) if (data2(i,2)>=0 && data2(i,2)<=20) table2(2,1)=table2(2,1)+1; elseif (data2(i,2)>=21 && data2(i,2)<=40) table2(2,2)=table2(2,2)+1; elseif (data2(i,2)>=41 && data2(i,2)<=60) table2(2,3)=table2(2,3)+1; elseif (data2(i,2)>=61 && data2(i,2)<=80) table2(2,4)=table2(2,4)+1; elseif (data2(i,2)>=81 && data2(i,2)<=100) table2(2,5)=table2(2,5)+1; end end end %%total for i=1:9 table(10,i)=sum(table(1:9,i)); end for i=1:9 table(i,10)=sum(table(i,1:9)); end
for i=1:5 table2(3,i)=sum(table2(1:2,i)); end for i=1:2 table2(i,6)=sum(table2(i,1:5)); end Step4: Calculate the marginal distributions, joint distributions, conditional distributions
Calculate the marginal, joint, conditional distributions according to their definitions.
%% marginal distribution
for i=1:9
pf(i)=table(10,i)/196;
end
for i=1:9
px(i)=table(i,10)/196;
end
for i=1:5
pf2(i)=table2(3,i)/196;
end
for i=1:2
px2(i)=table2(i,6)/196;
end
%% joint distribution
for i=1:9
for j= 1:9
p_joint(i,j)=table(i,j)/196;
end
end
for i=1:2
for j= 1:5
p_joint2(i,j)=table2(i,j)/196;
end
end
%conditional distribution p(f|x)
for i=1:9
for j= 1:9
p_fx(i,j)=p_joint(i,j)/px(i);
end
end
for i=1:2
for j= 1:5
p_fx2(i,j)=p_joint2(i,j)/px2(i);
end
end
%conditional distribution p(x|f)
for i=1:9
for j= 1:9
p_xf(i,j)=p_joint(i,j)/pf(j);
end
end
for i=1:2
for j= 1:5
p_xf2(i,j)=p_joint2(i,j)/pf2(j);
end
end Step 5: PlottingDon't forget to create a new figure (i.e. figure1=figure) for each plot plot marginal distributions
 %marginal distribution for temperature
figure1=figure;
h1=subplot(2,1,1,'Parent',figure1,...
'XTick',[1:1:9],'XTickLabel',[10:10:100],...
'FontSize',10,'FontWeight','bold');
hold(h1,'all');
plot(px(:),'-.','LineWidth',2,'Color',[0 0.8 0.8]);
xlim([1 9]);
title('observation','FontSize',12,'FontWeight','bold');
xlabel('temperature(F)','FontSize',8,'FontWeight','bold');
ylabel('marginal probability','FontSize',8,'FontWeight','bold');
set(gca, ...
'Box' , 'on' , ...
'TickDir' , 'in' , ...
'TickLength' , [.01 .02] , ...
'XMinorTick' , 'off' , ...
'YMinorTick' , 'off' , ...
'YGrid' , 'off' , ...
'XColor' , [.3 .3 .3], ...
'YColor' , [.3 .3 .3],...
'LineWidth' , 2 );
h2=subplot(2,1,2,'Parent',figure1,...
'XTick',[1:1:9],'XTickLabel',[10:10:100],...
'FontSize',10,'FontWeight','bold');
hold(h2,'all');
plot(pf(:),'-.','LineWidth',2,'Color',[0 0.8 0.8]);
xlim([1 9]);
title('forecasts','FontSize',12,'FontWeight','bold');
xlabel('temperature(F)','FontSize',8,'FontWeight','bold');
ylabel('marginal probability','FontSize',8,'FontWeight','bold');
set(gca, ...
'Box' , 'on' , ...
'TickDir' , 'in' , ...
'TickLength' , [.01 .02] , ...
'XMinorTick' , 'off' , ...
'YMinorTick' , 'off' , ...
'YGrid' , 'off' , ...
'XColor' , [.3 .3 .3], ...
'YColor' , [.3 .3 .3],...
'LineWidth' , 2 );
%marginal distribution for precipitation
figure2=figure;
h1=subplot(2,1,1,'Parent',figure2,...
'XTick',[1:1:2],'XTickLabel',[0:100:100],...
'FontSize',10,'FontWeight','bold');
hold(h1,'all');
plot(px2(:),'-.','LineWidth',2,'Color',[0 0.8 0.8]);
xlim([1 2]);
title('observation','FontSize',12,'FontWeight','bold');
xlabel('POP','FontSize',8,'FontWeight','bold');
ylabel('marginal probability','FontSize',8,'FontWeight','bold');
set(gca, ...
'Box' , 'on' , ...
'TickDir' , 'in' , ...
'TickLength' , [.01 .02] , ...
'XMinorTick' , 'off' , ...
'YMinorTick' , 'off' , ...
'YGrid' , 'off' , ...
'XColor' , [.3 .3 .3], ...
'YColor' , [.3 .3 .3],...
'LineWidth' , 2 );
h2=subplot(2,1,2,'Parent',figure2,...
'XTick',[1:1:5],'XTickLabel',[0:20:100],...
'FontSize',10,'FontWeight','bold');
hold(h2,'all');
plot(pf2(:),'-.','LineWidth',2,'Color',[0 0.8 0.8]);
xlim([1 5]);
title('forecasts','FontSize',12,'FontWeight','bold');
xlabel('POP','FontSize',8,'FontWeight','bold');
ylabel('marginal probability','FontSize',8,'FontWeight','bold');
set(gca, ...
'Box' , 'on' , ...
'TickDir' , 'in' , ...
'TickLength' , [.01 .02] , ...
'XMinorTick' , 'off' , ...
'YMinorTick' , 'off' , ...
'YGrid' , 'off' , ...
'XColor' , [.3 .3 .3], ...
'YColor' , [.3 .3 .3],...
'LineWidth' , 2 ); Plot the scatter points plots
%scatter points plot for temperature
figure3=figure;
scatter(data(:,1),data(:,2));
xlim([0 100]);
ylim([0 100]);
xlabel('observations','FontWeight','bold','FontSize',12);
ylabel('forecasts','FontWeight','bold','FontSize',12);
set(gca, ...
'Box' , 'on' , ...
'TickDir' , 'in' , ...
'TickLength' , [.01 .02] , ...
'XMinorTick' , 'off' , ...
'YMinorTick' , 'off' , ...
'YGrid' , 'off' , ...
'XColor' , [.3 .3 .3], ...
'YColor' , [.3 .3 .3],...
'LineWidth' , 2 );
%scatter points plot for precipitation
figure4=figure;
scatter(data2(:,1),data2(:,2));
xlim([-50 150]);
ylim([-50 150]);
xlabel('observations','FontWeight','bold','FontSize',12);
ylabel('forecasts','FontWeight','bold','FontSize',12);
set(gca, ...
'Box' , 'on' , ...
'TickDir' , 'in' , ...
'TickLength' , [.01 .02] , ...
'XMinorTick' , 'off' , ...
'YMinorTick' , 'off' , ...
'YGrid' , 'off' , ...
'XColor' , [.3 .3 .3], ...
'YColor' , [.3 .3 .3],...
'LineWidth' , 2 ); Plot absolute error as a function of f
%plot the absolute error as a function of f, when p is set to differnt
%value
figure5=figure;
for i = 1:101
x(i)= (i-1)*0.01;
p=0;
y0(i)=p+(1-2*p)*x(i);
p=0.2;
y2(i)=p+(1-2*p)*x(i);
p=0.4;
y4(i)=p+(1-2*p)*x(i);
p=0.5;
y5(i)=p+(1-2*p)*x(i);
p=0.6;
y6(i)=p+(1-2*p)*x(i);
p=0.8;
y8(i)=p+(1-2*p)*x(i);
p=1.0;
y10(i)=p+(1-2*p)*x(i);
end
plot(x,y0,'LineWidth',2);
hold all;
plot(x,y2,'LineWidth',2);
plot(x,y4,'LineWidth',2);
plot(x,y5,'LineWidth',2);
plot(x,y6,'LineWidth',2);
plot(x,y8,'LineWidth',2);
plot(x,y10,'LineWidth',2);
hleg=legend('p=0','p=0.2','p=0.4','p=0.5',...
'p=0.6','p=0.8','p=1.0',...
'Orientation','horizontal');
set(hleg,'TextColor',[.3,.3,.3],'FontSize',8);
xlim([0 1]);
ylim([0 1]);
xlabel('forecast','FontWeight','bold','FontSize',12);
ylabel('E[AE]','FontWeight','bold','FontSize',12);
set(gca, ...
'Box' , 'on' , ...
'TickDir' , 'in' , ...
'TickLength' , [.01 .02] , ...
'XMinorTick' , 'off' , ...
'YMinorTick' , 'off' , ...
'YGrid' , 'off' , ...
'XColor' , [.3 .3 .3], ...
'YColor' , [.3 .3 .3],...
'LineWidth' , 2 );
|