Search in sources :

Example 1 with Variate

use of dr.stats.Variate in project beast-mcmc by beast-dev.

the class SkylineReconstructor method main.

public static void main(String[] argv) {
    Variate x = null;
    List<Variate> plots = new ArrayList<Variate>();
    for (int i = 1; i <= 200; i++) {
        String stem = "sim" + (i < 10 ? "00" : (i < 100 ? "0" : "")) + i;
        try {
            SkylineReconstructor skyline = new SkylineReconstructor(new File(stem + ".log"), new File(stem + ".trees"), 1000000, 200, 0.0, 150000, 0.0);
            if (x == null) {
                x = skyline.getXData();
            }
            plots.add(skyline.getYDataMean());
        } catch (IOException e) {
            e.printStackTrace();
        } catch (ImportException e) {
            e.printStackTrace();
        } catch (TraceException e) {
            e.printStackTrace();
        }
        if (i % 10 == 0) {
            System.err.println("Read " + i);
        }
    }
    for (int i = 0; i < x.getCount(); i++) {
        System.out.print(x.get(i));
        for (Variate y : plots) {
            System.out.print("\t" + y.get(i));
        }
        System.out.println();
    }
}
Also used : ImportException(jebl.evolution.io.ImportException) TraceException(dr.inference.trace.TraceException) Variate(dr.stats.Variate) ArrayList(java.util.ArrayList) IOException(java.io.IOException) File(java.io.File)

Example 2 with Variate

use of dr.stats.Variate in project beast-mcmc by beast-dev.

the class LineageCountThroughTime method getLTT.

public static Variate[] getLTT(String treeFile, double minTime, double maxTime, int binCount, // number of trees to skip
int skip) throws IOException, ImportException {
    double delta = (maxTime - minTime) / (binCount - 1);
    BufferedReader reader = new BufferedReader(new FileReader(treeFile));
    String line = reader.readLine();
    TreeImporter importer;
    if (line.toUpperCase().startsWith("#NEXUS")) {
        importer = new NexusImporter(reader);
    } else {
        importer = new NewickImporter(reader, false);
    }
    int state = 0;
    while (importer.hasTree() && state < skip) {
        importer.importNextTree();
        state += 1;
    }
    // the age of the end of this group
    List<double[]> branchingTimes = new ArrayList<double[]>();
    state = 0;
    while (importer.hasTree()) {
        RootedTree tree = (RootedTree) importer.importNextTree();
        double[] bt = new double[tree.getInternalNodes().size()];
        int i = 0;
        for (Node node : tree.getInternalNodes()) {
            bt[i] = tree.getHeight(node);
            i++;
        }
        Arrays.sort(bt);
        branchingTimes.add(bt);
        state += 1;
    }
    Variate[] bins = new Variate[binCount];
    double height = 0;
    double n = branchingTimes.get(0).length;
    for (int k = 0; k < binCount; k++) {
        bins[k] = new Variate.D();
        if (height >= 0.0 && height <= maxTime) {
            for (state = 0; state < branchingTimes.size(); state++) {
                int index = 0;
                while (index < branchingTimes.get(state).length && branchingTimes.get(state)[index] < height) {
                    index += 1;
                }
                double lineageCount = 1;
                if (index < branchingTimes.get(state).length) {
                    lineageCount = n - index + 1;
                }
                bins[k].add(lineageCount);
            }
        }
        height += delta;
    }
    Variate xData = new Variate.D();
    Variate yDataMean = new Variate.D();
    Variate yDataMedian = new Variate.D();
    Variate yDataUpper = new Variate.D();
    Variate yDataLower = new Variate.D();
    double t = minTime;
    for (Variate bin : bins) {
        xData.add(t);
        if (bin.getCount() > 0) {
            yDataMean.add(bin.getMean());
            yDataMedian.add(bin.getQuantile(0.5));
            yDataLower.add(bin.getQuantile(0.025));
            yDataUpper.add(bin.getQuantile(0.975));
        } else {
            yDataMean.add(Double.NaN);
            yDataMedian.add(Double.NaN);
            yDataLower.add(Double.NaN);
            yDataUpper.add(Double.NaN);
        }
        t += delta;
    }
    return new Variate[] { xData, yDataMean };
}
Also used : NexusImporter(jebl.evolution.io.NexusImporter) Variate(dr.stats.Variate) Node(jebl.evolution.graphs.Node) ArrayList(java.util.ArrayList) RootedTree(jebl.evolution.trees.RootedTree) NewickImporter(jebl.evolution.io.NewickImporter) BufferedReader(java.io.BufferedReader) TreeImporter(jebl.evolution.io.TreeImporter) FileReader(java.io.FileReader)

Example 3 with Variate

use of dr.stats.Variate in project beast-mcmc by beast-dev.

the class KernelDensityEstimate method calculateDensity.

protected void calculateDensity(Variate data, int minimumBinCount) {
    FrequencyDistribution frequency = calculateFrequencies(data, minimumBinCount);
    xCoordinates = new Variate.D();
    yCoordinates = new Variate.D();
    double x = frequency.getLowerBound() - (frequency.getBinSize() / 2.0);
    int extraEdgeCount = 0;
    while (kde.pdf(x) > minDensity && x > lowerBoundary) {
        x -= frequency.getBinSize();
        extraEdgeCount += 1;
    }
    xCoordinates.add(x);
    yCoordinates.add(0.0);
    x += frequency.getBinSize();
    int count = 0;
    while (count < (frequency.getBinCount() + extraEdgeCount)) {
        // ||
        //                (kde.pdf(x) > minDensity && x < upperBoundary)) {
        xCoordinates.add(x);
        yCoordinates.add(kde.pdf(x));
        x += frequency.getBinSize();
        count++;
    }
    //        System.err.println("kde = " + kde.pdf(x));
    while (kde.pdf(x) > minDensity) {
        //            System.err.println("add bit on end!!!");
        xCoordinates.add(x);
        yCoordinates.add(kde.pdf(x));
        x += frequency.getBinSize();
    }
    xCoordinates.add(x);
    yCoordinates.add(0.0);
//
//
//        int extraBinsOnEdges = 5;
//        double x = frequency.getLowerBound() - extraBinsOnEdges * frequency.getBinSize();
//        for (int i = 0; i < frequency.getBinCount() + 2 * extraBinsOnEdges; i++) {
//            double xMidPoint = x + (frequency.getBinSize() / 2.0);
//            xData.add(xMidPoint);
//            yData.add(kde.pdf(xMidPoint));
//            x += frequency.getBinSize();
//        }
}
Also used : Variate(dr.stats.Variate) FrequencyDistribution(dr.util.FrequencyDistribution)

Aggregations

Variate (dr.stats.Variate)3 ArrayList (java.util.ArrayList)2 TraceException (dr.inference.trace.TraceException)1 FrequencyDistribution (dr.util.FrequencyDistribution)1 BufferedReader (java.io.BufferedReader)1 File (java.io.File)1 FileReader (java.io.FileReader)1 IOException (java.io.IOException)1 Node (jebl.evolution.graphs.Node)1 ImportException (jebl.evolution.io.ImportException)1 NewickImporter (jebl.evolution.io.NewickImporter)1 NexusImporter (jebl.evolution.io.NexusImporter)1 TreeImporter (jebl.evolution.io.TreeImporter)1 RootedTree (jebl.evolution.trees.RootedTree)1