%$RCSfile: DNASequenzEval.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)


disp(' ... weiter mit Taste...');pause
disp('Teilaufgabe (b)');
LL=[10:10:90 100:100:1000]; results=zeros(length(LL),5);
for l=1:length(LL)
  DNAcounts=CountACGT(DNASequenz(1:LL(l)));
  results(l,:) = [LL(l), ...
		  MultinomialDistribution(M25,  DNAcounts), ...
		  MultinomialDistribution(M50,  DNAcounts), ...
		  MultinomialDistribution(M100, DNAcounts), ...
		  MultinomialDistribution(M0,   DNAcounts)];
end
results(:,2:5) = log10(results(:,2:5));
disp(results); plot(results(:,1), results(:,2:5))
%gset term postscript;gset output "DNASequenzEvalB.ps";replot;gset term x11

disp(' ... weiter mit Taste...');pause
disp('Teilaufgabe (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,:));
end
Lp2=length(p2);
disp([p2(1:5:Lp2); Likelihood(:,1:5:Lp2)]')
plot(p2, Likelihood)
%gset term postscript;gset output "DNASequenzEvalC.ps";replot;gset term x11
