830 - MarkovModel::Impute()

void MarkovModel::Impute(char * major, char * observed, float * probs,
                         char ** haplotypes, float ** freqs, int position)
   {
   double P[5] = {0.0, 0.0, 0.0, 0.0, 0.0};

   for (int i = 0; i < states; i++)
      P[haplotypes[i][position]] += probs[i];

   double ptotal = P[1] + P[2] + P[3] + P[4];
   double pmajor = P[major[position]];

   int mle = 1;
   for (int i = 2; i <= 4; i++)
      if (P[i] >= P[mle])
         mle = i;

   char labels[] = {0, 'a', 'c', 'g', 't'};

   imputedDose[position] += imputedHap[position] = (pmajor / ptotal);
   imputedAlleles[position] = labels[mle];

   // What would we have imputed if we didn't observe a genotype at the current position?
   // This assumes that freqs[i][position] stores the frequency of allele i at position (for i > 0)
   // and is equal to 1.0 when the current genotype is not observed (i = 0)

   double fmatch = 1.0 / (1. - E[position] + E[position] * freqs[observed[position]][position]);
   double fmismatch = 1.0 / (E[position] * freqs[observed[position]][position]);

   for (int i = 1; i <= 4; i++)
      if (observed[position] == i)
         P[i] *= fmatch;
      else
         P[i] *= fmismatch;

   ptotal = P[1] + P[2] + P[3] + P[4];
   pmajor = P[major[position]];

   leaveOneOut[position] = pmajor / ptotal;
   }