White-noise clock model for topology estimation

[davidcerny]

Lepage et al. (2007) described the white-noise clock model, in which each branch of the tree draws its rate multiplier from a gamma distribution with a mean of 1 and a variance that is inversely proportional to that branch's temporal duration. As a result, rates of temporally long branches tend to stick close to the global mean rate, while the rates of temporally short branches can oscillate around that mean rate wildly.

For many years, this model was available in MrBayes under the name "IGR" ("independent gamma rates"), before the developers concluded that this was somewhat misleading, as it implies a model where individual branch rates are independent draws from a single gamma distribution (referred to as "UGam" by Lepage et al. 2007), similar to the well-known uncorrelated lognormal model (UCLN). Accordingly, in the current development version of MrBayes ("v3.2.8"), IGR refers to the i.i.d. model (UGam), while WN refers to what IGR used to mean in the previous versions.

I came up with a RevBayes implementation of the white-noise model based on the UGam code shown in the "Dating with relative constraints" tutorial. You can see it in the mvFossilTimeSlideUniform_multiple_runs.Rev test script, which is in turn based on a script I used in Pol et al. (2024):

for (i in n_branches:1) {
    # Get branch duration in calendar time units
    times[i] = fbd_tree.branchLength(i)

    # Relative branch-wise rates: these serve as multipliers to the global_rate_mean (c).
    # Each relative rate r_i is gamma distributed with mean = 1.0 and variance equal
    # to sigma/(t_i * c). Since dnGamma() is parameterized using shape and rate in
    # RevBayes, this corresponds to a shape of (t_i * c)/sigma and a rate equal to the
    # shape.
        
    gamma_param[i] := (times[i] * global_rate_mean)/sigma
    rel_branch_rates[i] ~ dnGamma(gamma_param[i], gamma_param[i])
    moves.append( mvScale(rel_branch_rates[i], lambda=0.5, weight=35/n_branches,
                          tune=true) )

    # The actual branch rate is the product of the global_rate_mean and the relative rate.
    branch_rates[i] := global_rate_mean * rel_branch_rates[i]
}

However, after I gave this code another look in light of #958 and #963, it occurred to me that this, too, was likely an unsuccessful attempt to convert code intended to be run on a fixed-topology to something that can be used for topology estimation. Since fbd_tree is a stochastic variable when the tree is estimated, we should be using the deterministic assignment operator := rather than = when extracting its branch lengths. In contrast to the tricky case seen in #958, this doesn't introduce a cycle into the model graph, nor does it require us to index by deterministic nodes (which is currently impossible).

So my question is – is this change (i.e., times[i] := fbd_tree.branchLength(i)) enough to turn the code above into a valid implementation of the white-noise model in RevBayes, or are there other problems with it that I'm overlooking?