From: Fritz Heinrichmeyer <fritz.heinrichmeyer@...>
Date: 2009-04-09T16:25:07+09:00
Subject: Re: what is happening here

define lambda, high and low as floats, lambda = 1.0, high=2.0, low = 0.0


Am Mittwoch, den 08.04.2009, 07:57 -0500 schrieb George George:
> > It seems like you're using the ** operator wrong.  You should use Ruby's
> > Math.exp instead (which is the direct equivalent to Perl's exp).
> > 
> > -Matthias
> 
> Thank you Mat
> i have changed that. here is the full script in ruby:
> Pn = 0.25 #probabilty of any nuclotide
> 
> match = 10
> mismatch = -10
> 
> expected_score = match * 0.25 + mismatch * 0.75
>  if match <=0 or expected_score >= 0
>    puts "illegal scores"
> #   quit
>  end
> 
> #calculate lambda
> #initial estimates
> lambda = 1
> high      = 2
> low        = 0
> 
> 
> while (high - low > 0.001) do
>   #calculate the sum of all normalized scores
>   sum = Pn * Pn * Math.exp(lambda * match) * 4 + Pn * Pn * 
> Math.exp(lambda * mismatch) * 12
> 
>   #refine guess at lambda
>   if sum > 1
>     high = lambda
>     lambda = (lambda  + low)/2
>   else
>     low = lambda
>     lambda = (lambda+high)/2
>   end
> end
> 
>   #compute target frequencies and H
>   target_id = Pn * Pn * Math.exp(lambda * match) * 4
>    h = lambda * match * target_id + lambda * mismatch * (1 - target_id)
> 
> puts "expected_score #{expected_score}\n"
> puts "lambda: #{lambda} nats, #{lambda/Math.log(2)} bits"
> puts "H: #{h} nats , #{h/Math.log(2)} bits"
> puts "%ID: #{target_id * 100}"
> 
> 
> and the original perl script is here:
> #!/usr/bin/perl -w use strict;
> use constant Pn => 0.25; # probability of any nucleotide
> 
> die "usage: $0 <match> <mismatch>\n" unless @ARGV == 2;
> 
> my ($match, $mismatch) = @ARGV;
> 
> my $expected_score = $match * 0.25 + $mismatch * 0.75;
> die "illegal scores\n" if $match <= 0 or $expected_score >= 0;
> 
> 
> # calculate lambda
> my ($lambda, $high, $low) = (1, 2, 0); # initial estimates
> 
> while ($high - $low > 0.001) {         # precision
>     # calculate the sum of all normalized scores
>     my $sum = Pn * Pn * exp($lambda * $match)    * 4
>             + Pn * Pn * exp($lambda * $mismatch) * 12;
>     # refine guess at lambda
>     if ($sum > 1) {
>         $high = $lambda;
>         $lambda = ($lambda + $low)/2;
>     }
>     else {
>         $low = $lambda;
>         $lambda = ($lambda + $high)/2;
>     }
> }
> 
> 
> # compute target frequency and H
> my $targetID = Pn * Pn * exp($lambda * $match) * 4;
> my $H = $lambda * $match    *     $targetID
>       + $lambda * $mismatch * (1 -$targetID);
> 
> 
> # output
> print "expscore: $expected_score\n";
> print "lambda:   $lambda nats (", $lambda/log(2), " bits)\n";
> print "H:        $H nats (", $H/log(2), " bits)\n";
> print "%ID:      ", $targetID * 100, "\n";
> 
> 
> The ruby version seems to get into an infinite loop. while the perl 
> version executes nicely and with the expected results.
>