Qing Zhu-Matlab how to

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 data
Open 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 tables
Two 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: Plotting
Don'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         );



ċ
data.txt
(4k)
Michael Baldwin,
Sep 9, 2012, 4:48 PM
v.4
Comments