I am implementing a categorical frequency matrix in HYPHY using the ConstructCategoryMatrix command.

However I am having some trouble interpreting the resulting matrix. Do the values  for a particular codon still need to be normalized?

Thanks!

Hi there,

ConstructCategoryMatrix takes a lot of options which control what it generates; one of the options will output normalized posterior probabilities. I'll post an example batch file later today.

Sergei

Thank you for the quick reply.

I will check out your batch file when you put it up. However, I have one quick question to make sure I understand the output from the ConstructCategoryMatrix command.

For a CategoryMatrix of {2,1}, if it gives you {0.1} {0.3}, then there is a distribution of 25% and 75% between the two categories.  

Thanks again.

Hi there,

By default, the function returns conditional likelihoods: Pr (site | rate class). In your example, the probability of observing a site given the model in class 1 is 0.1, and given the model in class 2 is 0.3

Most of the time one wants Pr (rate class | site), which can be obtained by a simple application of the Bayes rule:


Code (]
Pr (rate class | site) = Pr (site|rate class) * Pr (rate class) / Pr (site)
[/code):

Pr (site) is simply the sum of all Pr (site|rate class) * Pr (rate class) (the normalization factor), and Pr (rate class) can be obtained by calling [code] ConstructCategoryMatrix (priors, lfID, WEIGHTS);

HTH,
Sergei

PS I could probably give you more pointers if I had a better idea of what the final objective was.

Hi Sergei,

The following code functions correctly, though I had to use for loops instead of matrix functions.

Is there a way to do it faster?

Thanks for the help,
John


Code (]ConstructCategoryMatrix(catmat,theLnLik,COMPLETE);
ConstructCategoryMatrix (priors, theLnLik, WEIGHTS);

result = {Rows(catmat),Columns(catmat)};

for( i=0; i<Rows(catmat); i=i+1)
{
 for(j=0; j<Columns(catmat); j=j+1)
 {
   result[i):

[j] = catmat[i][j]*priors[i][0];  } } for( j=0; j<Columns(catmat); j=j+1) {  sum = 0;  for( i=0; i<Rows(catmat); i=i+1)  {    sum = sum+result[i][j];  }  for( i=0; i<Rows(catmat); i=i+1)  {    result[i][j] = result[i][j]/sum;  } } fprintf(stdout,"\nresult:",result,"\n");

Hi John,

There really is not a built-in way to do it much faster (in terms of lines of code). HyPhy template batch files contain some HBL functions to do the same, but your code is perfectly adequate. You can define it as a function and include it in other files for later use.

Sergei