Welcome, Guest. Please Login
YaBB - Yet another Bulletin Board
 
  HomeHelpSearchLogin  
 
setting constraints (Read 2439 times)
Rachel
YaBB Newbies
*
Offline


Feed your monkey!

Posts: 9
setting constraints
Mar 19th, 2009 at 3:06pm
 
Hello,

I’ve been taking a closer look at the ins and outs of setting constraints and the output from my test files is not what I expect.  

I set up things as follows:


DataSet myData = ReadDataFile("/Users/rachel/HYPHY/test.seq");

DataSetFilter ElementFilter = CreateFilter(myData, 1, "0-71");
DataSetFilter ReferenceFilter = CreateFilter(myData, 1, "72-11366");

HarvestFrequencies(obsFreqsRef, ReferenceFilter, 1, 1, 1);
HarvestFrequencies(obsFreqsEle, ElementFilter, 1, 1, 1);


FMatrix = {{*,mu,mu,mu}{mu,*,mu,mu}{mu,mu,*,mu}{mu,mu,mu,*}};

Model FRef = (FMatrix, obsFreqsRef);
Tree myTreeRef = (((hg18, panTro2), rheMac2), (mm8, rn4), canFam2);

Model FEle = (FMatrix, obsFreqsEle);
Tree myTreeEle = (((hg18, panTro2), rheMac2), (mm8, rn4), canFam2);

I build and optimize the likelihood function:

LikelihoodFunction theLnLik = (ReferenceFilter, myTreeRef, ElementFilter, myTreeEle);
Optimize (params, theLnLik);


I now apply the constraints.  

global G=0;
ReplicateConstraint(“this1.?.mu:=G*this2.?.mu”, myTreeEle, myTreeRef);
Optimize (params, theLnLik);

I get (in part) the following output:

G=0.57508
myTreeRef.hg18.mu=0.00597841
myTreeEle.hg18.mu=G*myTreeRef.hg18.mu=0.00343806
myTreeEle.Node2.mu=G*myTreeRef.Node2.mu=0.00980094

I then add another constraint:

ClearConstraints(myTreeEle.hg18.mu);
global S:>1;
myTreeEle.hg18.mu := S*G*myTreeRef.hg18.mu;
Optimize (params, theLnLik);
:

I get (in part) the following output:

G=0.482008
S=32.4987
myTreeRef.hg18.mu=0.00548192
myTreeEle.hg18.mu=S*G*myTreeRef.hg18.mu=0.0858723
myTreeEle.Node2.mu=G*myTreeRef.Node2.mu=0.00822833

It seems that G and the branch lengths have been re-optimized.   I would like to hold G and the branch lengths constant (at the values obtained under the first set of constraints) and only optimize S.  Is there a way to do this and am I correct in my assessment of what is going on?  


Thanks very much in advance for your help.

Cheers,

Rachel
Back to top
 
 
IP Logged
 
Sergei
YaBB Administrator
*****
Offline


Datamonkeys are forever...

Posts: 1658
UCSD
Gender: male
Re: setting constraints
Reply #1 - Mar 19th, 2009 at 6:37pm
 
Dear Rachel,

Can you attach the batch file to the message because some of the characters were lost in copy/paste.

Cheers,
Sergei
Back to top
 

Associate Professor
Division of Infectious Diseases
Division of Biomedical Informatics
School of Medicine
University of California San Diego
WWW WWW  
IP Logged
 
Rachel
YaBB Newbies
*
Offline


Feed your monkey!

Posts: 9
Re: setting constraints
Reply #2 - Mar 20th, 2009 at 10:13am
 
Hi Sergei,

Sorry about that.  It looks ok on my browser. 

The batch file is attached to this message. 

Thanks

Rachel
Back to top
 
Multimedia File Viewing and Clickable Links are available for Registered Members only!!  You need to Login Login (1 KB | )
 
IP Logged
 
Sergei
YaBB Administrator
*****
Offline


Datamonkeys are forever...

Posts: 1658
UCSD
Gender: male
Re: setting constraints
Reply #3 - Mar 20th, 2009 at 5:11pm
 
Dear Rachel,

After

Code:
ClearConstraints(myTreeEle.hg18.mu);
 



add

Code:
G:=G__; /*constraint G to the current value; __ is the in-place evaluation construct */
ReplicateConstraint ("this1.?.?:=this2.?.?__", myTreeRef, myTreeRef);
/* constrain branch lengths in myTreeRef to current values */
 



Cheers,
Sergei
Back to top
 

Associate Professor
Division of Infectious Diseases
Division of Biomedical Informatics
School of Medicine
University of California San Diego
WWW WWW  
IP Logged
 
Rachel
YaBB Newbies
*
Offline


Feed your monkey!

Posts: 9
Re: setting constraints
Reply #4 - Mar 20th, 2009 at 5:36pm
 
That answers my question.  Thanks very much!

Rachel
Back to top
 
 
IP Logged