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]);
}
}
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;
}
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();
}
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);
}
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());
}
Aggregations