perl/etc/justonce/20140324hw3.pl

   1 #!/usr/bin/perl -w
   2 use strict;
   3 use warnings;
   4 use Data::Dump;
   5
   6 die "Usage: $0 <s> <t>\nW(AA), W(Aa) and W(aa) is 1-s, 1, 1-t respectively.\n" if @ARGV<2;
   7
   8 my ($s,$t)=@ARGV;
   9
  10 my ($w11,$w12,$w22) = (1-$s, 1, 1-$t);
  11 my ($AA,$Aa,$aa) = (0.50, 0.40, 0.10);
  12
  13 my $p = $AA + $Aa/2;
  14 my $q = $aa + $Aa/2;
  15
  16 my $avgW = $p*$p*$w11 + 2*$p*$q*$w12 + $q*$q*$w22;
  17
  18 print "GT: ($AA,$Aa,$aa) => ($p,$q)\nFit: ($s,$t) => ($w11,$w12,$w22), $avgW\n";
  19
  20 for my $i ( 1 .. 2 ) {
  21         my $avgW = $p*$p*$w11 + 2*$p*$q*$w12 + $q*$q*$w22;
  22         my $np = ( $p*$p*$w11 + $p*$q*$w12 ) / $avgW;
  23         my $nq = ( $p*$q*$w12 + $q*$q*$w22 ) / $avgW;
  24         $AA = $p*$p; $Aa = 2*$p*$q; $aa = $q*$q;
  25         print "$i\t-> ($np, $nq) ",$np+$nq,"\t($AA,$Aa,$aa)\tavg(w)=$avgW\n";
  26         $p = $np; $q = $nq;
  27 }