Hard-coded exp priors on IOU parameters

src/init.c

@@ -3579,7 +3579,7 @@ void PHYREX_Set_Default_Migrep_Mod(int n_otu, t_phyrex_mod *t)
   t->model_id        = -1;
   t->use_locations   = -1;
   t->sampling_scheme = -1;
-  t->safe_phyrex     = NO;
+  t->safe_phyrex     = YES;
   t->dist_type       = HAVERSINE;
 
   t->lim_up->lonlat[0] = 100.;
@@ -3604,9 +3604,6 @@ void PHYREX_Set_Default_Migrep_Mod(int n_otu, t_phyrex_mod *t)
   t->prior_param_rad   = 1.;
   t->update_rad        = NO;
 
-  t->min_sigsq         = 0.0;
-  t->max_sigsq         = 1.E+3;
-  t->prior_param_sigsq = 10.0;
 
   t->min_veloc         = -1.E+5;
   t->max_veloc         = +1.E+5;
@@ -3615,16 +3612,19 @@ void PHYREX_Set_Default_Migrep_Mod(int n_otu, t_phyrex_mod *t)
   t->min_omega         = 0.0;
   t->max_omega         = 1.E+3;
 
-  t->ou_theta             = 1.0;
-  t->min_ou_theta         = 1.0E-4;
-  t->max_ou_theta         = 10.0;
+  t->ou_theta             = 1.0E-1;
+  t->min_ou_theta         = 0.0;
+  t->max_ou_theta         = 1000.0;
 
   for(int i=0;i<t->n_dim;++i) t->ou_mu[i] = 0.0;
   t->min_ou_mu            = -10.;
   t->max_ou_mu            = 10.;
 
+  t->min_sigsq         = 0.0;
+  t->max_sigsq         = 1.E+3;
+  t->prior_param_sigsq = 10.0;
   assert(t->n_dim > 0);
-  for(int i=0;i<t->n_dim;++i) t->sigsq[i] = t->min_sigsq + (t->max_sigsq-t->min_sigsq)/20.;
+  for(int i=0;i<t->n_dim;++i) t->sigsq[i] = 1.E-1;
 
   t->nu = 1.0E-0;
 

src/iou.c

@@ -65,7 +65,8 @@ phydbl IOU_Velocity_Variance_Along_Edge(t_node *d, short int dim, t_tree *tree)
     {
       PhyML_Printf("\n. ta: %f t: %f sg2: %f sinh: %f",ta,t,sg2,sinh(ta*t));
     }
-
+//   PhyML_Printf("\n. t: %f ta: %f sg2: %f var: %f",t,ta,sg2,var);
+
   assert(isnan(var) == NO);
   assert(isinf(var) == NO);
 
@@ -115,22 +116,24 @@ phydbl IOU_Location_Mean_Along_Edge(t_node *d, short int dim, t_tree *tree)
    at both extremities */
 phydbl IOU_Location_Variance_Along_Edge(t_node *d, short int dim, t_tree *tree)
 {
-  phydbl var,t,ta,sg;
+  phydbl var,t,ta,sg2;
 
   assert(d != tree->n_root);
 
 
   t  = fabs(tree->times->nd_t[d->num] - tree->times->nd_t[d->anc->num]);
   ta = tree->mmod->ou_theta;
-  sg = tree->mmod->sigsq[dim];
+  sg2 = tree->mmod->sigsq[dim] * tree->mmod->sigsq_scale[d->num] * tree->mmod->sigsq_scale_norm_fact;
 
   if(t < SMALL) return(0.0);
 
-  /* var = sg/pow(ta,3)*(ta*t-2.*(cosh(ta*t)-1.)/sinh(ta*t)); */
-  var = sg/pow(ta,3)*(ta*t-2.*(1./tanh(ta*t) - 1./sinh(ta*t)));
+//   var = sg/pow(ta,3)*(ta*t-2.*(cosh(ta*t)-1.)/sinh(ta*t));
+  var = sg2/pow(ta,3)*(ta*t-2.*(1./tanh(ta*t) - 1./sinh(ta*t)));
+
+  if(var > sg2*pow(t,3)/12) var = sg2*pow(t,3)/12;
 
   if(var < 0.0) var = 0.0;
-
+    
   assert(isnan(var) == NO);
 
   return(var);
@@ -156,9 +159,9 @@ phydbl IOU_Prior_Theta(t_tree *tree)
 {
   phydbl lbda,lnP;
 
-  lbda = 1.0;
-  lnP = log(lbda)-lbda*tree->mmod->ou_theta;
-
+  lnP = 0.0;
+  lbda = 1.0E+4;
+  lnP += Log_Dexp(tree->mmod->ou_theta,lbda);
   return(lnP);
 }
 
@@ -167,7 +170,13 @@ phydbl IOU_Prior_Theta(t_tree *tree)
 
 phydbl IOU_Prior_Mu(t_tree *tree)
 {
-  return(0.0);
+  phydbl lbda,lnP;
+
+  lnP = 0.0;
+  lbda = 1.0E+3;
+  for (int i = 0; i < tree->mmod->n_dim; i++)
+      lnP += Log_Dexp(tree->mmod->ou_mu[i], lbda);
+  return (lnP);
 }
 
 //////////////////////////////////////////////////////////////