[beast-mcmc] r6740 committed - Gibbs-disguised operators to sample from gamma and normal distribution...

[beast...@googlecode.com]

Revision: 6740
Author: bael...@gmail.com
Date: Thu Apr 16 14:37:42 2015 UTC
Log: Gibbs-disguised operators to sample from gamma and normal
distributions.
https://code.google.com/p/beast-mcmc/source/detail?r=6740

Added:

/trunk/src/dr/inference/operators/GibbsIndependentJointNormalGammaOperator.java
Modified:
/trunk/src/dr/inference/operators/GibbsIndependentGammaOperator.java

=======================================
--- /dev/null
+++
/trunk/src/dr/inference/operators/GibbsIndependentJointNormalGammaOperator.java
Thu Apr 16 14:37:42 2015 UTC
@@ -0,0 +1,238 @@
+/*
+ * GibbsIndependentNormalDistributionOperator.java
+ *
+ * Copyright (c) 2002-2014 Alexei Drummond, Andrew Rambaut and Marc Suchard
+ *
+ * This file is part of BEAST.
+ * See the NOTICE file distributed with this work for additional
+ * information regarding copyright ownership and licensing.
+ *
+ * BEAST is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * BEAST is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with BEAST; if not, write to the
+ * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
+ * Boston, MA 02110-1301 USA
+ */
+
+package dr.inference.operators;
+
+import cern.jet.random.Gamma;
+import cern.jet.random.Normal;
+import cern.jet.random.engine.MersenneTwister;
+import cern.jet.random.engine.RandomEngine;
+import dr.inference.distribution.NormalDistributionModel;
+import dr.inference.model.Bounds;
+import dr.inference.model.Parameter;
+import dr.inference.model.Variable;
+import dr.math.MathUtils;
+import dr.math.distributions.GammaDistribution;
+import dr.xml.*;
+
+/**
+ * An independent joint normal gamma sampler to propose new (independent)
values from provided normal and gamma distribution models.
+ *
+ * @author Guy Baele
+ *
+ */
+public class GibbsIndependentJointNormalGammaOperator extends
SimpleMCMCOperator implements GibbsOperator {
+
+ public static final String OPERATOR_NAME
= "GibbsIndependentJointNormalGammaOperator";
+ public static final String MEAN = "mean";
+ public static final String PREC = "precision";
+ public static final String SHAPE = "shape";
+ public static final String SCALE = "scale";
+
+ private Variable<Double> mean = null;
+ private Variable<Double> precision = null;
+ private NormalDistributionModel model = null;
+ private GammaDistribution gamma = null;
+ private boolean updateAllIndependently = true;
+
+ private static final boolean TRY_COLT = true;
+ private static final boolean DEBUG = false;
+ private static RandomEngine randomEngine;
+ private static Normal coltNormal;
+ //private static Gamma coltGamma;
+
+ public GibbsIndependentJointNormalGammaOperator(Variable mean, Variable
precision, NormalDistributionModel model, GammaDistribution gamma) {
+
+ this(mean, precision, model, gamma, 1.0);
+
+ }
+
+ public GibbsIndependentJointNormalGammaOperator(Variable mean, Variable
precision, NormalDistributionModel model, GammaDistribution gamma, double
weight) {
+
+ this(mean, precision, model, gamma, weight, true);
+
+ }
+
+ public GibbsIndependentJointNormalGammaOperator(Variable mean, Variable
precision, NormalDistributionModel model, GammaDistribution gamma, double
weight, boolean updateAllIndependently) {
+
+ this.mean = mean;
+ this.precision = precision;
+ this.model = model;
+ this.gamma = gamma;
+ this.updateAllIndependently = updateAllIndependently;
+ setWeight(weight);
+
+ if (TRY_COLT) {
+ randomEngine = new MersenneTwister(MathUtils.nextInt());
+ //create standard normal distribution, internal states will be
bypassed anyway
+ //takes mean and standard deviation
+ coltNormal = new Normal(0.0, 1.0, randomEngine);
+ //coltGamma = new Gamma(gamma.getShape(),
1.0/gamma.getScale(), randomEngine);
+ } else {
+ //no random draw with specified mean and stdev implemented in the
normal distribution in BEAST (as far as I know)
+ throw new RuntimeException("Normal distribution in BEAST still
needs a random sampler.");
+ }
+
+ }
+
+ public String getPerformanceSuggestion() {
+ return "";
+ }
+
+ public String getOperatorName() {
+ return "GibbsIndependentJointNormalGamma(" + mean.getVariableName()
+ "," + precision.getVariableName() + ")";
+ }
+
+ public int getStepCount() {
+ return 1;
+ }
+
+ /**
+ * change the parameter and return the hastings ratio.
+ */
+ public double doOperation() throws OperatorFailedException {
+
+ //double logq = 0;
+
+ //double currentValue;
+ double newValue;
+
+ final Bounds<Double> meanBounds = mean.getBounds();
+ final Bounds<Double> precBounds = precision.getBounds();
+ final int dim = mean.getSize();
+
+ if (updateAllIndependently) {
+ for (int i = 0; i < dim; i++) {
+
+ if (DEBUG) {
+ System.out.println("old precision value: " +
precision.getValue(i));
+ System.out.println("old mean value: " +
mean.getValue(i));
+ System.out.println("model mean check: " +
model.getMean().getValue(i));
+ System.out.println("model precision check: " +
model.getPrecision().getValue(i) + "\n");
+ }
+ //start with updating the precision
+ newValue = gamma.nextGamma();
+ //newValue = coltGamma.nextDouble();
+ while (newValue == 0.0) {
+ newValue = gamma.nextGamma();
+ //newValue = coltGamma.nextDouble();
+ }
+
+ if (newValue < precBounds.getLowerLimit(i) || newValue >
precBounds.getUpperLimit(i)) {
+ throw new OperatorFailedException("proposed value from
gamma distribution outside boundaries");
+ }
+
+ precision.setValue(i, newValue);
+
+ if (DEBUG) {
+ System.out.println("new precision value: " +
precision.getValue(i));
+ System.out.println("old mean value: " +
mean.getValue(i));
+ System.out.println("model mean check: " +
model.getMean().getValue(i));
+ System.out.println("model precision check: " +
model.getPrecision().getValue(i) + "\n");
+ }
+
+ //use the current mean and precision (standard deviation)
+ newValue = coltNormal.nextDouble(model.getMean().getValue(i), 1.0 /
Math.sqrt(model.getPrecision().getValue(i)));
+
+ //newValue = (double)model.nextRandom();
+
+ //System.out.print("normal distribution model: N(" +
model.getMean().getValue(i) + "," + model.getPrecision().getValue(i)
+ ") ");
+ //System.out.println(1.0 /
Math.sqrt(model.getPrecision().getValue(i)));
+ //System.out.println("current value: " + currentValue + " -- new
value: " + newValue);
+
+ //logq += (model.logPdf(currentValue) - model.logPdf(newValue));
+
+ if (newValue < meanBounds.getLowerLimit(i) || newValue >
meanBounds.getUpperLimit(i)) {
+ throw new OperatorFailedException("proposed value from
normal distribution outside boundaries");
+ }
+
+ mean.setValue(i, newValue);
+
+ if (DEBUG) {
+ System.out.println("new precision value: " +
precision.getValue(i));
+ System.out.println("new mean value: " +
mean.getValue(i));
+ System.out.println("model mean check: " +
model.getMean().getValue(i));
+ System.out.println("model precision check: " +
model.getPrecision().getValue(i) + "\n");
+ }
+
+ }
+ }
+
+ //return logq;
+ return 0;
+
+ }
+
+ public static dr.xml.XMLObjectParser PARSER = new
dr.xml.AbstractXMLObjectParser() {
+
+ public String getParserName() {
+ return OPERATOR_NAME;
+ }
+
+ public Object parseXMLObject(XMLObject xo) throws XMLParseException {
+
+ double weight = xo.getDoubleAttribute(WEIGHT);
+ double shape = xo.getDoubleAttribute(SHAPE);
+ double scale = xo.getDoubleAttribute(SCALE);
+
+ NormalDistributionModel model = (NormalDistributionModel)
xo.getChild(NormalDistributionModel.class);
+
+ XMLObject cxo = xo.getChild(MEAN);
+ Parameter mean = (Parameter) cxo.getChild(Parameter.class);
+
+ cxo = xo.getChild(PREC);
+ Parameter precision = (Parameter)
cxo.getChild(Parameter.class);
+
+ return new GibbsIndependentJointNormalGammaOperator(mean, precision,
model, new GammaDistribution(shape, scale), weight);
+ }
+
+
//************************************************************************
+ // AbstractXMLObjectParser implementation
+
//************************************************************************
+
+ public XMLSyntaxRule[] getSyntaxRules() {
+ return rules;
+ }
+
+ private final XMLSyntaxRule[] rules = {
+ AttributeRule.newDoubleRule(WEIGHT),
+ AttributeRule.newDoubleRule(SHAPE),
+ AttributeRule.newDoubleRule(SCALE),
+ new ElementRule(NormalDistributionModel.class),
+ new ElementRule(MEAN, new XMLSyntaxRule[]{new
ElementRule(Parameter.class)}),
+ new ElementRule(PREC, new XMLSyntaxRule[]{new
ElementRule(Parameter.class)})
+ };
+
+ public String getParserDescription() {
+ return "This element returns an independence sampler, disguised as a
Gibbs operator, from provided normal and gamma distributions.";
+ }
+
+ public Class getReturnType() {
+ return MCMCOperator.class;
+ }
+
+ };
+
+}
=======================================
--- /trunk/src/dr/inference/operators/GibbsIndependentGammaOperator.java
Wed Dec 3 10:21:42 2014 UTC
+++ /trunk/src/dr/inference/operators/GibbsIndependentGammaOperator.java
Thu Apr 16 14:37:42 2015 UTC
@@ -162,7 +162,7 @@
};

public String getParserDescription() {
- return "This element returns an independence sampler, disguised as a
Gibss operator, from a provided gamma prior.";
+ return "This element returns an independence sampler, disguised as a
Gibbs operator, from a provided gamma prior.";
}

public Class getReturnType() {