alex
YaBB Newbies
Offline
Posts: 8
Gender:
|
After making those changes (at least my interpretation of those changes), the branch lengths of the two trees are still different, and they are different from those of the reference tree. However, they appear to be different from the reference tree by only a scaling factor that is the same for each branch within a tree.
Are these scaling factors being estimated? If so, how can I fix them to 1.0? On the other hand if they are deterministic consequences of, for example, normalizing the expected rate of the component rate matrices, then how exactly are they calculated?
I have a separate question about how the likelihood is maximized. Is HyPhy taking advantage of the 'mixture' structure of my likelihood function as it searches for the MLE? It is possible to analytically find the maximum likelihood values of some of the parameters in the model conditional on the the data and the values of the other parameters.
Thanks, Alex
[code] VERBOSITY_LEVEL = 11;
DataSet spectrinData = ReadDataFile ("sim2.nex");
DataSetFilter filteredData = CreateFilter (spectrinData,1);
HarvestFrequencies (observedFreqs, filteredData, 1, 1, 1);
fprintf (stdout, observedFreqs, "\n");
/* stationary distribution in nucleotide frequencies */ global eqFreqA = 0.25; global eqFreqC = 0.25; global eqFreqG = 0.25; global eqFreqT := 1.0 - eqFreqA - eqFreqC - eqFreqG; eqFreqT :> 0;
/* another distribution for nucleotide frequencies */ global eqFreqA2 = 0.25; global eqFreqC2 = 0.25; global eqFreqG2 = 0.25; global eqFreqT2 := 1.0 - eqFreqA2 - eqFreqC2 - eqFreqG2; eqFreqT2 :> 0;
/* transition/transversion ratio */ global kappa = 1.0; kappa :> 0;
HKY85RateMatrix = {{*,a,kappa*a,a} {a,*,a,kappa*a} {kappa*a,a,*,a} {a,kappa*a,a,*}};
/* construct vectors from global variables */ estFreqs1 = { {eqFreqA, eqFreqC, eqFreqG, eqFreqT} };
estFreqs2 = { {eqFreqA2, eqFreqC2, eqFreqG2, eqFreqT2} }; /* assign models to trees */ Model m1 = (HKY85RateMatrix, estFreqs1);
ACCEPT_BRANCH_LENGTHS = 1; Tree givenTree = DATAFILE_TREE; branchNames = BranchName (givenTree, -1); /* vector populated with names for each branch in tree */ branchLengths = BranchLength (givenTree, -1); /* vector populated with branch lengths */ for (k = 0; k < Columns(branchNames)-1; k = k+1) { ExecuteCommands("givenTree."+branchNames[k]+".a:="+branchLengths[k]+";"); }
Model m2 = (HKY85RateMatrix, estFreqs2);
ACCEPT_BRANCH_LENGTHS = 1; Tree otherTree = DATAFILE_TREE; branchNames = BranchName (otherTree, -1); /* vector populated with names for each branch in tree */ branchLengths = BranchLength (otherTree, -1); /* vector populated with branch lengths */ for (k = 0; k < Columns(branchNames)-1; k = k+1) { ExecuteCommands("otherTree."+branchNames[k]+".a:="+branchLengths[k]+";"); }
global P = 0.5; /* mixing proportion */ P :> 0; P :< 1;
/* construct likelihood function using mixture of two rate distributions */ LikelihoodFunction lf = (filteredData, givenTree, filteredData, otherTree, "Log(SITE_LIKELIHOOD[0]*P+(1-P)*SITE_LIKELIHOOD[1])");
Optimize(res, lf);
fprintf(stdout, res); fprintf(stdout, lf);
[/code]
|