%$RCSfile: aufg_3_3.m,v $ $Revision: 1.1.1.1 $
format compact

DNASequenz; % lade DNASequenz
M0=1/4* [1 1 1 1]; % das richtige Modell

disp('Teilaufgabe (a)');
M25  = CountACGT(DNASequenz(1:25));  M25  = M25 /sum(M25)
M50  = CountACGT(DNASequenz(1:50));  M50  = M50 /sum(M50)
M100 = CountACGT(DNASequenz(1:100)); M100 = M100/sum(M100)
M500 = CountACGT(DNASequenz(1:500)); M500 = M500/sum(M500)
M1000= CountACGT(DNASequenz(1:1000));M1000= M1000/sum(M1000)

M= [M50;M100;M500;M1000];

figure(0);
disp(' ... weiter mit Taste...');pause
disp("\n");
disp('Teilaufgabe (b)');
disp("\nPlot zu p_a = p_t = 0.5, p_g = 0.5 - p_c...");
p1=1/4; p2=[0:0.005:0.5]; p3=0.5-p2; p4=1/4;
LL=[50 100 500 1000];Likelihood=zeros(length(LL), length(p2));
for iL=1:length(LL)
  DNAcounts=CountACGT(DNASequenz(1:LL(iL)));
  for i=1:length(p2)
    Likelihood(iL,i)=MultinomialDistribution([p1 p2(i) p3(i) p4], DNAcounts);
  end  
    Likelihood(iL,:)=Likelihood(iL,:)/max(Likelihood(iL,:));
  z = max(Likelihood(iL,:));
  y = find(z == Likelihood(iL,:));
  printf("N = %d, Maximum fuer p_c = %f\n", LL(iL), p2(y));
end
Lp2=length(p2);
#disp([p2(1:5:Lp2); Likelihood(:,1:5:Lp2)]')
title("p_a = p_t = 0.5");
plot(p2, Likelihood)
hold on;
x = [M(4,2),M(4,2)];
y = [0,1];
plot(x,y);
hold off;
gset term postscript;gset output "aufg_3_3b_falsch.ps";replot;gset term x11

figure(1);
disp("\nZum Vergleich: \"korrekter\" Plot, mit den geschaetzten Werten aus a)...");
LL=[50 100 500 1000];Likelihood=zeros(length(LL), length(p2));
for iL=1:length(LL)
  p1= M(iL,1); p4=M(iL,4);
  rest= 1 - p1 - p4;
  p2=[0:0.005:rest]; p3=rest-p2; 
  DNAcounts=CountACGT(DNASequenz(1:LL(iL)));
  for i=1:length(p2)
    Likelihood(iL,i)=MultinomialDistribution([p1 p2(i) p3(i) p4], DNAcounts);
  end  
    Likelihood(iL,:)=Likelihood(iL,:)/max(Likelihood(iL,:));
  z = max(Likelihood(iL,:));
  y = find(z == Likelihood(iL,:));
  printf("N = %d, Maximum fuer p_c = %f\n", LL(iL), p2(y));
end
Lp2=length(p2);
#disp([p2(1:5:Lp2); Likelihood(:,1:5:Lp2)]')
title("p_a und p_t aus Teil a)");
plot(p2, Likelihood)
hold on;
x = [M(4,2),M(4,2)];
y = [0,1];
plot(x,y);
hold off;
gset term postscript;gset output "aufg_3_3b_richtig.ps";replot;gset term x11
