%% Plot
x=[1:100]/100;
H1=log2((x).^(-x))+log2(((1-x)).^(-(1-x)));
H2=log2((x).^(-x))+2*log2((((1-x)/2)).^(-(1-x)/2));
H3=log2((x).^(-x))+3*log2((((1-x)/3)).^(-(1-x)/3));
H4=log2((x).^(-x))+4*log2((((1-x)/4)).^(-(1-x)/4));
%% initial 
show_ini= 1;
if show_ini
    tic;
    for letter = 1:nCONS+nVOWS  % can change the number of letters to create fewer charts
    %for letter = 1:7
    %     iSub = logical(eval(['i' ,regexprep(letters{letter},'(\<[a-z])','${upper($1)}')]));
        figure
        for i = 1:15
            if letter == 3
                Ninit(letter,i) = 1; % to account for letter G
            end
            if i > length(RD)
                plot(errorxsub_initial(letter,i),entropyxsub_initial(letter,i),'xb',...
                    'Marker',markers{16-i},'MarkerSize',100/sqrt(Ninit(letter,i)));hold on;
            else
                plot(errorxsub_initial(letter,i),entropyxsub_initial(letter,i),'xr',...
                    'Marker',markers{i},'MarkerSize',100/sqrt(Ninit(letter,i)));hold on;
            end;
        end
        l=legend('Alina','Angela','Edward','Latisha','Laura','Norene',...
            'Shauna','Teddy','Tony','Anton','Bob','Carly','Evan',...
            'Joanna','Miguel');
        set(l,'Location','EastOutside', 'FontSize',12);
        %legend('boxoff')
        %set(l,'Position',[0.5 0.15 0.035 0.05]);
        title(['H(p_e) for Letter (' letters{letter} ') in Initial Position'])
        ylabel('H [bits]')
        xlabel('p_e')
        % plot reference lines
%         plot(x,H1,'LineStyle','-','Color',[0.6 0.6 0.6]);hold on;
%         plot(1-x,H2,'LineStyle','-','Color',[0.6 0.6 0.6]);
%         plot(1-x,H3,'LineStyle','-','Color',[0.6 0.6 0.6]);
%         plot(1-x,H4,'LineStyle','-','Color',[0.6 0.6 0.6]);
        % label data points
        labels = SUBJECTS;
        ixSub = (1:15);
        
        pRD(letter,1:3) = polyfit(errorxsub_initial(letter,1:length(RD)),...
        entropyxsub_initial(letter,1:length(RD)),2)
        x1 = 0:0.001:(max(errorxsub_initial(letter,1:length(RD)))+.1);
        y1 = polyval(pRD(letter,1:3),x1);
        hold on
        plot(x1,y1,'b')
        pRC(letter,1:3) = polyfit(errorxsub_initial(letter,(length(RD)+1):end),...
            entropyxsub_initial(letter,(length(RD)+1):end),2)
        x2 = 0:0.001:(max(errorxsub_initial(letter,1:length(RD)))+.1);
        y2 = polyval(pRC(letter,1:3),x2);
        hold on
        plot(x2,y2,'r')
        if letter ~= 3
        axis([0,max(errorxsub_initial(letter,1:length(RD)))+.1, 0, ...
            2.5])%max(entropyxsub_initial(letter,1:length(RD)))+.25])
        end
    end
end