From: George George <george.githinji@...>
Date: 2009-04-08T21:57:19+09:00
Subject: Re: what is happening here


> 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.

-- 
Posted via http://www.ruby-forum.com/.