Examples with Likelihood dr.inference.model.Likelihood used on opensource projects

Example 26 with Likelihood

use of dr.inference.model.Likelihood in project beast-mcmc by beast-dev.

the class DirichletProcessOperator method doOperate.

//END: doOp
private void doOperate() throws MathException {
    // int index = 0;
    for (int index = 0; index < realizationCount; index++) {
        int[] occupancy = new int[uniqueRealizationCount];
        for (int i = 0; i < realizationCount; i++) {
            if (i != index) {
                int j = (int) categoriesParameter.getParameterValue(i);
                occupancy[j]++;
            }
        // END: i check
        }
        if (DEBUG) {
            System.out.println("N[-index]: ");
            dr.app.bss.Utils.printArray(occupancy);
        }
        Likelihood clusterLikelihood = (Likelihood) likelihood.getLikelihood(index);
        //			Likelihood clusterLikelihood = likelihood;
        int category = (int) categoriesParameter.getParameterValue(index);
        double value = uniqueParameters.getParameterValue(category);
        double[] clusterProbs = new double[uniqueRealizationCount];
        for (int i = 0; i < uniqueRealizationCount; i++) {
            double logprob = 0;
            if (occupancy[i] == 0) {
                // draw new
                // draw from base model, evaluate at likelihood
                double candidate = dpp.baseModel.nextRandom()[0];
                uniqueParameters.setParameterValue(category, candidate);
                double loglike = clusterLikelihood.getLogLikelihood();
                uniqueParameters.setParameterValue(category, value);
                logprob = Math.log((intensity) / (realizationCount - 1 + intensity)) + loglike;
            } else {
                // draw existing
                // likelihood for component x_index
                double candidate = dpp.getUniqueParameter(i).getParameterValue(0);
                uniqueParameters.setParameterValue(category, candidate);
                double loglike = clusterLikelihood.getLogLikelihood();
                uniqueParameters.setParameterValue(category, value);
                logprob = Math.log(occupancy[i]) / (realizationCount - 1 + intensity) + loglike;
            }
            // END: occupancy check
            clusterProbs[i] = logprob;
        }
        // END: i loop
        dr.app.bss.Utils.exponentiate(clusterProbs);
        if (DEBUG) {
            System.out.println("P(z[index] | z[-index]): ");
            dr.app.bss.Utils.printArray(clusterProbs);
        }
        // sample
        int sampledCluster = MathUtils.randomChoicePDF(clusterProbs);
        categoriesParameter.setParameterValue(index, sampledCluster);
        if (DEBUG) {
            System.out.println("sampled category: " + sampledCluster + "\n");
        }
    }
// END: index loop
}

Example 27 with Likelihood

use of dr.inference.model.Likelihood in project beast-mcmc by beast-dev.

the class ManyUniformGeoDistributionModelParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    List<GeoSpatialDistribution> distributions = new ArrayList<GeoSpatialDistribution>();
    List<Parameter> parameters = new ArrayList<Parameter>();
    List<Likelihood> oldLikelihood = new ArrayList<Likelihood>();
    for (int i = 0; i < xo.getChildCount(); ++i) {
        Object cxo = xo.getChild(i);
        MultivariateDistributionLikelihood ab = null;
        if (cxo instanceof MultivariateDistributionLikelihood) {
            ab = (MultivariateDistributionLikelihood) cxo;
        } else if (cxo instanceof CachedDistributionLikelihood) {
            ab = (MultivariateDistributionLikelihood) ((CachedDistributionLikelihood) cxo).getDistributionLikelihood();
        }
        if (ab != null) {
            parameters.add(ab.getDataParameter());
            MultivariateDistribution md = ab.getDistribution();
            oldLikelihood.add((Likelihood) cxo);
            if (md instanceof GeoSpatialDistribution) {
                distributions.add((GeoSpatialDistribution) md);
            } else {
                throw new XMLParseException("Unhandled distribution type in '" + xo.getId() + "'");
            }
        }
    }
    return new ManyUniformGeoDistributionModel(xo.getId(), parameters, distributions, oldLikelihood);
}

Example 28 with Likelihood

use of dr.inference.model.Likelihood in project beast-mcmc by beast-dev.

the class MCMCParser method parseXMLObject.

/**
     * @return an mcmc object based on the XML element it was passed.
     */
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    MCMC mcmc = new MCMC(xo.getAttribute(NAME, "mcmc1"));
    long chainLength = xo.getLongIntegerAttribute(CHAIN_LENGTH);
    boolean useCoercion = xo.getAttribute(COERCION, true);
    long coercionDelay = chainLength / 100;
    if (xo.hasAttribute(PRE_BURNIN)) {
        coercionDelay = xo.getIntegerAttribute(PRE_BURNIN);
    }
    coercionDelay = xo.getAttribute(COERCION_DELAY, coercionDelay);
    double temperature = xo.getAttribute(TEMPERATURE, 1.0);
    long fullEvaluationCount = xo.getAttribute(FULL_EVALUATION, 2000);
    double evaluationTestThreshold = MarkovChain.EVALUATION_TEST_THRESHOLD;
    if (System.getProperty("mcmc.evaluation.threshold") != null) {
        evaluationTestThreshold = Double.parseDouble(System.getProperty("mcmc.evaluation.threshold"));
    }
    evaluationTestThreshold = xo.getAttribute(EVALUATION_THRESHOLD, evaluationTestThreshold);
    int minOperatorCountForFullEvaluation = xo.getAttribute(MIN_OPS_EVALUATIONS, 1);
    MCMCOptions options = new MCMCOptions(chainLength, fullEvaluationCount, minOperatorCountForFullEvaluation, evaluationTestThreshold, useCoercion, coercionDelay, temperature);
    OperatorSchedule opsched = (OperatorSchedule) xo.getChild(OperatorSchedule.class);
    Likelihood likelihood = (Likelihood) xo.getChild(Likelihood.class);
    likelihood.setUsed();
    if (Boolean.valueOf(System.getProperty("show_warnings", "false"))) {
        // check that all models, parameters and likelihoods are being used
        for (Likelihood l : Likelihood.FULL_LIKELIHOOD_SET) {
            if (!l.isUsed()) {
                java.util.logging.Logger.getLogger("dr.inference").warning("Likelihood, " + l.getId() + ", of class " + l.getClass().getName() + " is not being handled by the MCMC.");
            }
        }
        for (Model m : Model.FULL_MODEL_SET) {
            if (!m.isUsed()) {
                java.util.logging.Logger.getLogger("dr.inference").warning("Model, " + m.getId() + ", of class " + m.getClass().getName() + " is not being handled by the MCMC.");
            }
        }
        for (Parameter p : Parameter.FULL_PARAMETER_SET) {
            if (!p.isUsed()) {
                java.util.logging.Logger.getLogger("dr.inference").warning("Parameter, " + p.getId() + ", of class " + p.getClass().getName() + " is not being handled by the MCMC.");
            }
        }
    }
    ArrayList<Logger> loggers = new ArrayList<Logger>();
    for (int i = 0; i < xo.getChildCount(); i++) {
        Object child = xo.getChild(i);
        if (child instanceof Logger) {
            loggers.add((Logger) child);
        }
    }
    mcmc.setShowOperatorAnalysis(true);
    if (xo.hasAttribute(OPERATOR_ANALYSIS)) {
        mcmc.setOperatorAnalysisFile(XMLParser.getLogFile(xo, OPERATOR_ANALYSIS));
    }
    Logger[] loggerArray = new Logger[loggers.size()];
    loggers.toArray(loggerArray);
    java.util.logging.Logger.getLogger("dr.inference").info("\nCreating the MCMC chain:" + "\n  chainLength=" + options.getChainLength() + "\n  autoOptimize=" + options.useCoercion() + (options.useCoercion() ? "\n  autoOptimize delayed for " + options.getCoercionDelay() + " steps" : "") + (options.getFullEvaluationCount() == 0 ? "\n  full evaluation test off" : ""));
    mcmc.init(options, likelihood, opsched, loggerArray);
    MarkovChain mc = mcmc.getMarkovChain();
    double initialScore = mc.getCurrentScore();
    if (initialScore == Double.NEGATIVE_INFINITY) {
        String message = "The initial posterior is zero";
        if (likelihood instanceof CompoundLikelihood) {
            message += ": " + ((CompoundLikelihood) likelihood).getDiagnosis(2);
        } else {
            message += "!";
        }
        throw new IllegalArgumentException(message);
    }
    if (!xo.getAttribute(SPAWN, true))
        mcmc.setSpawnable(false);
    return mcmc;
}

Example 29 with Likelihood

use of dr.inference.model.Likelihood in project beast-mcmc by beast-dev.

the class CompoundGaussianProcessParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    List<GaussianProcessRandomGenerator> gpList = new ArrayList<GaussianProcessRandomGenerator>();
    List<Likelihood> likelihoodList = new ArrayList<Likelihood>();
    List<Integer> copyList = new ArrayList<Integer>();
    for (int i = 0; i < xo.getChildCount(); ++i) {
        Object obj = xo.getChild(i);
        GaussianProcessRandomGenerator gp = null;
        Likelihood likelihood = null;
        int copies = -1;
        if (obj instanceof DistributionLikelihood) {
            DistributionLikelihood dl = (DistributionLikelihood) obj;
            if (!(dl.getDistribution() instanceof GaussianProcessRandomGenerator)) {
                throw new XMLParseException("Not a Gaussian process");
            }
            likelihood = dl;
            gp = (GaussianProcessRandomGenerator) dl.getDistribution();
            copies = 0;
            for (Attribute<double[]> datum : dl.getDataList()) {
                //                    Double draw = (Double) gp.nextRandom();
                //                    System.err.println("DL: " + datum.getAttributeName() + " " + datum.getAttributeValue().length + " " + "1");
                copies += datum.getAttributeValue().length;
            }
        } else if (obj instanceof MultivariateDistributionLikelihood) {
            MultivariateDistributionLikelihood mdl = (MultivariateDistributionLikelihood) obj;
            if (!(mdl.getDistribution() instanceof GaussianProcessRandomGenerator)) {
                throw new XMLParseException("Not a Gaussian process");
            }
            likelihood = mdl;
            gp = (GaussianProcessRandomGenerator) mdl.getDistribution();
            copies = 0;
            double[] draw = (double[]) gp.nextRandom();
            for (Attribute<double[]> datum : mdl.getDataList()) {
                //                    System.err.println("ML: " + datum.getAttributeName() + " " + datum.getAttributeValue().length + " " + draw.length);
                copies += datum.getAttributeValue().length / draw.length;
            }
        } else if (obj instanceof GaussianProcessRandomGenerator) {
            gp = (GaussianProcessRandomGenerator) obj;
            likelihood = gp.getLikelihood();
            copies = 1;
        } else {
            throw new XMLParseException("Not a Gaussian process");
        }
        gpList.add(gp);
        likelihoodList.add(likelihood);
        copyList.add(copies);
    }
    //        System.exit(-1);
    return new CompoundGaussianProcess(gpList, likelihoodList, copyList);
}

Example 30 with Likelihood

use of dr.inference.model.Likelihood in project beast-mcmc by beast-dev.

the class BirthDeathSSLikelihoodTest method testBirthDeathLikelihoodBEAST2.

public void testBirthDeathLikelihoodBEAST2() {
    System.out.println("RootHeight = " + tree2.getRootHeight());
    Variable<Double> origin = new Variable.D("origin", 6.0);
    final double birthRate = 2.0;
    final double deathRate = 1.0;
    // rate of sampling taxa through time
    final double psiRate = 0.5;
    // the proportion of taxa sampled, default to fix to 0
    final double sampleProbability = 0.0;
    final boolean hasFinalSample = false;
    Variable<Double> b = new Variable.D("b", birthRate);
    Variable<Double> d = new Variable.D("d", deathRate);
    Variable<Double> psi = new Variable.D("psi", psiRate);
    Variable<Double> p = new Variable.D("p", sampleProbability);
    // sampleBecomesNonInfectiousProb
    Variable<Double> r = new Variable.D("r", 0.0);
    SpeciationModel speciationModel = new BirthDeathSerialSamplingModel(b, d, psi, p, false, r, hasFinalSample, origin, Units.Type.YEARS);
    Likelihood likelihood = new SpeciationLikelihood(tree2, speciationModel, "bdss.like");
    assertEquals(-19.0198, likelihood.getLogLikelihood(), 1e-5);
}