4ment / phylostan Goto Github PK

I am using the following commands to generate the stan script, and run the stan script respectively:

phylostan build -s sample.stan -m JC69 -C 4 --estimate_rate --clock strict --coalescent constant
phylostan run -s sample.stan -m JC69 -C 4 --estimate_rate --clock strict --coalescent constant -i out.fasta -t sample_newick.newick -o sample -q meanfield

I would expect that there are branch lengths from the bottommost tree in sample.trees, as given below, but, there are none.

I've also tried looking in the file populated in the -o flag of the run command, labelled as VB samples below. I have fewer than 10 sequences, yet there are 10 nonzero branch lengths in the VB samples file, so I am not sure what to make of this. Is there another place that I can get inferred branch lengths? Or is there something I'm doing incorrectly?

Original data
VB samples
Script generated using build
Tree topology
sample.trees

According to the Stan documentation, defining a density using _log (eg. real constant_coalescent_log(real[] heights, real popSize, int[,] map, real[] lowers){) is deprecated. That definition should become

real constant_coalescent_lpdf(real[] heights | real popSize, int[,] map, real[] lowers){

instead.

To do topology comparisons, we want to extract ELBOs from stan runs.
In cmdstan these are saved to a 'diagnostic file' that can be specified on the command-line, but pystan & rstan do not provide an interface to specify this.
There are two possible paths -- either wrap cmdstan from the calling python script, or improve the pystan API so that the diagnostic file can be given.
The former appears complicated.
The latter also seemed initially complicated, but investigation into the pystan codebase highlights a path. In particular, code needs to be added to open the diagnostic file around this line:
https://github.com/stan-dev/pystan/blob/develop/pystan/stan_fit.hpp#L1030
and 'diagnostic_file' needs to be listed as a valid argument around this line:
https://github.com/stan-dev/pystan/blob/develop/pystan/model.py#L868

there may be a few other little edits required -- but looks pretty straightforward and possibly something that could be sent as a pull request to the standev guys.

Hello,

Congrats on the excellent work. More than issues, mine are questions on how to use the software:

• I have not fully understood how to try phylostan for variable population size. I have tried in the two stages (build and run) to add "-c skyride" or "-c skygrid," but I always get somehow estimates for the constant population. At least, that's what the output on the terminal suggests. The output file also includes the same number of parameters, and there is a single theta (I would expect for variable population size this number to change)
• Could you also please share some post-processing files/instructions? I am unsure how to manage the output file, so I have not been able to generate any of the posterior plots myself.

Thanks a lot!

Lorenzo