view release on metacpan or  search on metacpan

bin/alnpi  view on Meta::CPAN

313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
## COMPUTE PROBABILITY OF ALLELIC CONFIGURATION USING KARLIN MACGREGOR 1972
$Num_alleles = scalar keys %alleles;
$H = 1;
my %a;
map {$H -= (($_/$nseqs) ** 2)} values %alleles;
if (!$slide_window and !$pairwise and $Theta_w > 0) {
  $Exp_num_alleles = 1 / $Theta_w;
  for my $i (1 .. ($#seqs )) {
    $Exp_num_alleles += 1 / ($Theta_w + $i);
  }
  $Exp_num_alleles *= $Theta_w; ## Ewen's Sampling formula
  map {$a{$_}++} values %alleles; ## THIS MAKES THE a_i
  my $log_prob   = $Num_alleles * log ($Theta_w);
  #      $log_prob  += Math::Gsl::Sf::sf_lnfact($Num_alleles);
  $log_prob  += Ramanujan_logfact($Num_alleles);
   
  foreach my $i  (0 .. ($nseqs - 1)) {
    $log_prob -= log ($Theta_w + $i);
  }
  foreach (keys %a) {
    #   $log_prob -= ($a{$_} * log ($_)) + Math::Gsl::Sf::sf_lnfact($a{$_});