Search in sources :

Example 16 with Vector

use of dr.math.matrixAlgebra.Vector in project beast-mcmc by beast-dev.

the class TwoStateMarkovRewardsTest method testTwoStateRewards.

public void testTwoStateRewards() {
    DataType dataType = TwoStates.INSTANCE;
    FrequencyModel freqModel = new FrequencyModel(TwoStates.INSTANCE, new double[] { 0.5, 0.5 });
    Parameter rates = new Parameter.Default(new double[] { 4.0, 6.0 });
    ComplexSubstitutionModel twoStateModel = new ComplexSubstitutionModel("two", dataType, freqModel, rates) {
    };
    twoStateModel.setNormalization(false);
    MarkovJumpsSubstitutionModel markovRewards = new MarkovJumpsSubstitutionModel(twoStateModel, MarkovJumpsType.REWARDS);
    double[] r = new double[2];
    double[] q = new double[4];
    double[] c = new double[4];
    int mark = 0;
    double weight = 1.0;
    r[mark] = weight;
    markovRewards.setRegistration(r);
    twoStateModel.getInfinitesimalMatrix(q);
    System.out.println("Q = " + new Vector(q));
    System.out.println("Reward for state 0");
    double time = 1.0;
    markovRewards.computeCondStatMarkovJumps(time, c);
    System.out.println("Reward conditional on X(0) = i, X(t) = j: " + new Vector(c));
    double endTime = 10.0;
    int steps = 10;
    for (time = 0.0; time < endTime; time += (endTime / steps)) {
        markovRewards.computeCondStatMarkovJumps(time, c);
        // start = 0, end = 0
        System.out.println(time + "," + c[0]);
    }
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) ComplexSubstitutionModel(dr.evomodel.substmodel.ComplexSubstitutionModel) DataType(dr.evolution.datatype.DataType) Parameter(dr.inference.model.Parameter) Vector(dr.math.matrixAlgebra.Vector) MarkovJumpsSubstitutionModel(dr.evomodel.substmodel.MarkovJumpsSubstitutionModel)

Example 17 with Vector

use of dr.math.matrixAlgebra.Vector in project beast-mcmc by beast-dev.

the class SumDerivative method getGradientLogDensity.

@Override
public double[] getGradientLogDensity() {
    int size = derivativeList.size();
    if (DEBUG) {
    // start timer
    }
    final double[] derivative = derivativeList.get(0).getGradientLogDensity();
    if (DEBUG) {
        String name = derivativeList.get(0).getLikelihood().getId();
        System.err.println(name);
        System.err.println(new Vector(derivative));
    }
    for (int i = 1; i < size; i++) {
        if (DEBUG) {
        // start timer
        }
        final double[] temp = derivativeList.get(i).getGradientLogDensity();
        if (DEBUG) {
            String name = derivativeList.get(i).getLikelihood().getId();
            System.err.println(name);
            System.err.println(new Vector(temp));
        }
        for (int j = 0; j < temp.length; j++) {
            derivative[j] += temp[j];
        }
    }
    if (DEBUG) {
        System.exit(-1);
    }
    return derivative;
}
Also used : Vector(dr.math.matrixAlgebra.Vector)

Example 18 with Vector

use of dr.math.matrixAlgebra.Vector in project beast-mcmc by beast-dev.

the class HiddenMarkovRatesTest method testHiddenRates.

public void testHiddenRates() {
    final int index = 0;
    double[] freqs = getBinaryFreqs(index);
    FrequencyModel binaryFreqModel = new FrequencyModel(TwoStates.INSTANCE, freqs);
    int relativeTo = 0;
    Parameter ratesParameter = new Parameter.Default(0);
    GeneralSubstitutionModel binaryModel = new GeneralSubstitutionModel("binary", TwoStates.INSTANCE, binaryFreqModel, ratesParameter, relativeTo);
    UniformizedSubstitutionModel uSM = new UniformizedSubstitutionModel(binaryModel, MarkovJumpsType.REWARDS);
    uSM.setSaveCompleteHistory(true);
    double[] rewardRegister = new double[] { 0.0, 1.0 };
    uSM.setRegistration(rewardRegister);
    final double[] hkyFreqs = getHKYFreqs(index);
    FrequencyModel hkyFreqModel = new FrequencyModel(Nucleotides.INSTANCE, hkyFreqs);
    final double kappa = getKappa(index);
    final HKY hky = new HKY(kappa, hkyFreqModel);
    final double length = getLength(index);
    double[] resultCompleteHistory = new double[16];
    final int replicates = getNumberReplicates(index);
    double result = 0.0;
    for (int r = 0; r < replicates; ++r) {
        result += oneCompleteHistoryReplicate(resultCompleteHistory, hky, uSM, length);
    }
    result /= replicates;
    normalize(resultCompleteHistory, replicates);
    System.out.println("Averaged probabilities");
    System.out.println(result);
    System.out.println(new Vector(resultCompleteHistory));
    System.out.println();
    double[] intermediate = new double[16];
    hky.getTransitionProbabilities(result, intermediate);
    System.out.println("Intermediate using above average reward");
    System.out.println(result);
    System.out.println(new Vector(intermediate));
    System.out.println();
    double[] resultExpected = new double[16];
    UniformizedSubstitutionModel expectedUSM = new UniformizedSubstitutionModel(binaryModel, MarkovJumpsType.REWARDS, replicates);
    expectedUSM.setRegistration(rewardRegister);
    result = oneCompleteHistoryReplicate(resultExpected, hky, expectedUSM, length);
    System.out.println("Averaged reward");
    System.out.println(result);
    System.out.println(new Vector(resultExpected));
    System.out.println();
    double[] originalProbs = new double[16];
    hky.getTransitionProbabilities(length, originalProbs);
    System.out.println("Original probabilities");
    System.out.println(new Vector(originalProbs));
    System.out.println();
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) HKY(dr.evomodel.substmodel.nucleotide.HKY) Parameter(dr.inference.model.Parameter) GeneralSubstitutionModel(dr.evomodel.substmodel.GeneralSubstitutionModel) Vector(dr.math.matrixAlgebra.Vector) UniformizedSubstitutionModel(dr.evomodel.substmodel.UniformizedSubstitutionModel)

Example 19 with Vector

use of dr.math.matrixAlgebra.Vector in project beast-mcmc by beast-dev.

the class TN93Test method testTN93.

public void testTN93() {
    Parameter kappa1 = new Parameter.Default(5.0);
    Parameter kappa2 = new Parameter.Default(2.0);
    double[] pi = new double[] { 0.40, 0.20, 0.30, 0.10 };
    double time = 0.1;
    FrequencyModel freqModel = new FrequencyModel(Nucleotides.INSTANCE, pi);
    TN93 tn = new TN93(kappa1, kappa2, freqModel);
    EigenDecomposition decomp = tn.getEigenDecomposition();
    Vector eval = new Vector(decomp.getEigenValues());
    System.out.println("Eval = " + eval);
    double[] probs = new double[16];
    tn.getTransitionProbabilities(time, probs);
    System.out.println("new probs = " + new Vector(probs));
    // check against old implementation
    dr.oldevomodel.substmodel.FrequencyModel oldFreq = new dr.oldevomodel.substmodel.FrequencyModel(Nucleotides.INSTANCE, pi);
    dr.oldevomodel.substmodel.TN93 oldTN = new dr.oldevomodel.substmodel.TN93(kappa1, kappa2, oldFreq);
    double[] oldProbs = new double[16];
    oldTN.getTransitionProbabilities(time, oldProbs);
    System.out.println("old probs = " + new Vector(oldProbs));
    assertEquals(probs, oldProbs, 10E-6);
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) TN93(dr.evomodel.substmodel.nucleotide.TN93) EigenDecomposition(dr.evomodel.substmodel.EigenDecomposition) Parameter(dr.inference.model.Parameter) Vector(dr.math.matrixAlgebra.Vector)

Example 20 with Vector

use of dr.math.matrixAlgebra.Vector in project beast-mcmc by beast-dev.

the class MarkovModulatedGY94CodonModel method main.

public static void main(String[] args) {
    GY94CodonModel codonModel = new GY94CodonModel(Codons.UNIVERSAL, new Parameter.Default(1.0), new Parameter.Default(2.0), new FrequencyModel(Codons.UNIVERSAL, new Parameter.Default(61, 1.0 / 61.0)));
    EigenDecomposition ed1 = codonModel.getEigenDecomposition();
    //        double[][] q = codonModel.getQ();
    //        System.err.println("matrixQ = \n"+codonModel.printQ());// new Matrix(q));
    FrequencyModel freqModel = new FrequencyModel(HiddenCodons.UNIVERSAL_HIDDEN_2, new Parameter.Default(122, 1.0 / 122.0));
    System.err.println("freq = " + new Vector(freqModel.getFrequencies()));
    //        System.exit(-1);
    MarkovModulatedGY94CodonModel mmCodonModel = new MarkovModulatedGY94CodonModel(HiddenCodons.UNIVERSAL_HIDDEN_2, new Parameter.Default(2, 5.0), new Parameter.Default(2, 1.0), new Parameter.Default(2.0), freqModel);
    EigenDecomposition ed2 = mmCodonModel.getEigenDecomposition();
    // new Matrix(q));
    System.err.println("matrixQ = \n" + mmCodonModel.printQ());
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) EigenDecomposition(dr.evomodel.substmodel.EigenDecomposition) Parameter(dr.inference.model.Parameter) Vector(dr.math.matrixAlgebra.Vector)

Aggregations

Vector (dr.math.matrixAlgebra.Vector)39 Matrix (dr.math.matrixAlgebra.Matrix)12 FrequencyModel (dr.evomodel.substmodel.FrequencyModel)9 Parameter (dr.inference.model.Parameter)6 NodeRef (dr.evolution.tree.NodeRef)5 WishartSufficientStatistics (dr.math.distributions.WishartSufficientStatistics)5 MarkovJumpsSubstitutionModel (dr.evomodel.substmodel.MarkovJumpsSubstitutionModel)4 HKY (dr.evomodel.substmodel.nucleotide.HKY)4 EigenDecomposition (dr.evomodel.substmodel.EigenDecomposition)3 DataType (dr.evolution.datatype.DataType)2 StateHistory (dr.inference.markovjumps.StateHistory)2 MultivariateNormalDistribution (dr.math.distributions.MultivariateNormalDistribution)2 IllegalDimension (dr.math.matrixAlgebra.IllegalDimension)2 CompleteHistorySimulator (dr.app.beagle.tools.CompleteHistorySimulator)1 ComplexSubstitutionModel (dr.evomodel.substmodel.ComplexSubstitutionModel)1 GeneralSubstitutionModel (dr.evomodel.substmodel.GeneralSubstitutionModel)1 ProductChainFrequencyModel (dr.evomodel.substmodel.ProductChainFrequencyModel)1 UniformizedSubstitutionModel (dr.evomodel.substmodel.UniformizedSubstitutionModel)1 TN93 (dr.evomodel.substmodel.nucleotide.TN93)1 ContinuousDiffusionIntegrator (dr.evomodel.treedatalikelihood.continuous.cdi.ContinuousDiffusionIntegrator)1