i have some perl code 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; } } trying to rewrite this in ruby while (high - low > 0.001) #calculate the sum of all normalized scores sum = Pn * Pn **(lambda * match) * 4 + Pn * Pn **(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 Is this the right way? i seem to get different results when i execute the ruby version. any ideas?

on 2009-04-08 14:01

on 2009-04-08 14:32

George George wrote: > sum = Pn * Pn **(lambda * match) * 4 + Pn * Pn **(lambda * mismatch) * > 12 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

on 2009-04-08 14:57

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

on 2009-04-08 15:25

in ruby lambda, high and low are integers and divide by two is also integer arithmetic. Perl uses float for scalars by default. Am Mittwoch, den 08.04.2009, 07:57 -0500 schrieb George George:

on 2009-04-08 15:32

Quoting George George <george.githinji@gmail.com>: > } > else { > $low = $lambda; > $lambda = ($lambda + $high)/2; > } > } > Pn * exp($lambda * $match) is Pn * (Math::E ** (lambda + match)), not Pn ** (lambda + match). WBR, Peter Zotov

on 2009-04-09 09:27

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:

on 2009-04-09 10:33

```
Fritz Heinrichmeyer wrote:
> define lambda, high and low as floats, lambda = 1.0, high=2.0, low = 0.0
Thank you so much everyone for your responses, it works nicely now!
GG
```