Skip to content

Commit

Permalink
fix Cauchy gradient and add junit
Browse files Browse the repository at this point in the history
  • Loading branch information
msuchard committed Oct 5, 2024
1 parent ea377bc commit 538e916
Show file tree
Hide file tree
Showing 3 changed files with 84 additions and 2 deletions.
1 change: 1 addition & 0 deletions src/dr/app/beast/development_parsers.properties
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@ dr.inferencexml.distribution.RowDimensionPoissonPriorParser
dr.inferencexml.distribution.MomentDistributionModelParser
dr.inferencexml.distribution.DeterminentalPointProcessPriorParser
dr.inferencexml.distribution.TruncatedDistributionLikelihoodParser
dr.inferencexml.distribution.CauchyDistributionModelParser

# STRUCTURED COALESCENT
dr.evomodel.coalescent.structure.StructuredCoalescentLikelihood
Expand Down
4 changes: 2 additions & 2 deletions src/dr/inference/distribution/CauchyDistribution.java
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ public CauchyDistribution(double median, double scale){

public static double pdf(double x, double median, double scale) {
double dev = x - median;
return scale / (dev * dev + scale * scale) / Math.PI;
return scale / (dev * dev + scale * scale) / Math.PI;
}

public static double cdf(double x, double median, double scale) {
Expand All @@ -69,7 +69,7 @@ public static double quantile(double p, double median, double scale) {

public static double gradLogPdf(double x, double median, double scale) {
double dev = x - median;
return 2 * dev / (dev * dev + scale * scale);
return -2 * dev / (dev * dev + scale * scale);
}

public static double logPdf(double x, double median, double scale) {
Expand Down
81 changes: 81 additions & 0 deletions src/test/dr/inference/distribution/CauchyDistributionTest.java
Original file line number Diff line number Diff line change
@@ -0,0 +1,81 @@
/*
* CauchyDistributionTest.java
*
* Copyright © 2002-2024 the BEAST Development Team
* http://beast.community/about
*
* 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 test.dr.inference.distribution;

import dr.inference.distribution.CauchyDistributionModel;
import dr.inference.model.Parameter;
import dr.math.NumericalDerivative;
import dr.math.UnivariateFunction;
import junit.framework.TestCase;
import org.apache.commons.math.distribution.CauchyDistributionImpl;

/**
* @author Marc A Suchard
*/

public class CauchyDistributionTest extends TestCase {

public CauchyDistributionTest(String name) {
super(name);
}

public void setUp() throws Exception {
super.setUp();
}

public void testCauchyDistribution() {
CauchyDistributionModel d1 = new CauchyDistributionModel(
new Parameter.Default(2.0),
new Parameter.Default(3.0)
);

CauchyDistributionImpl d2 = new CauchyDistributionImpl(2, 3);

assertEquals(d1.pdf(-1), d2.density(-1), 1E-10);
assertEquals(d1.logPdf(-1), Math.log(d2.density(-1)), 1E-10);

double dx = NumericalDerivative.firstDerivative(new UnivariateFunction() {
@Override
public double evaluate(double argument) {
return d1.logPdf(argument);
}

@Override
public double getLowerBound() {
return Double.NEGATIVE_INFINITY;
}

@Override
public double getUpperBound() {
return Double.POSITIVE_INFINITY;
}
}, 3.0);

assertEquals(d1.getGradientLogDensity(new double[] { 3.0 } )[0], dx, 1E-4);
}
}

0 comments on commit 538e916

Please sign in to comment.