impl BBMRF getQ()

beast-dev · Nov 4, 2023 · 54fe8b5 · 54fe8b5
1 parent b980090
commit 54fe8b5
Show file tree

Hide file tree

Showing 2 changed files with 22 additions and 5 deletions.
diff --git a/src/dr/math/distributions/BayesianBridgeMarkovRandomField.java b/src/dr/math/distributions/BayesianBridgeMarkovRandomField.java
@@ -55,10 +55,27 @@ public BayesianBridgeMarkovRandomField(String name,
     }
 
     @Override
-    SymmetricTriDiagonalMatrix getQ() {
+    protected SymmetricTriDiagonalMatrix getQ() {
         if (!qKnown) {
-//            bayesianBridge.getStandardDeviation(i);
-            throw new RuntimeException("Not yet implemented");
+
+            final double[] diagonal = Q.diagonal;
+            final double[] offDiagonal = Q.offDiagonal;
+
+            double sd = bayesianBridge.getStandardDeviation(0);
+            offDiagonal[0] = -1.0 / (sd * sd);
+            diagonal[0] = -offDiagonal[0];
+
+            for (int i = 1; i < dim - 1; ++i) {
+                sd = bayesianBridge.getStandardDeviation(i);
+                offDiagonal[i] = -1.0 / (sd * sd);
+                diagonal[i] = -(offDiagonal[i - 1] + offDiagonal[i]);
+            }
+
+            diagonal[dim - 1] = -offDiagonal[dim - 2];
+
+            // TODO Update for lambda != 1 and for weights
+
+            qKnown = true;
         }
         return Q;
     }

diff --git a/src/dr/math/distributions/GaussianMarkovRandomField.java b/src/dr/math/distributions/GaussianMarkovRandomField.java
@@ -43,7 +43,7 @@ public class GaussianMarkovRandomField extends RandomFieldDistribution {
 
     public static final String TYPE = "GaussianMarkovRandomField";
 
-    private final int dim;
+    protected final int dim;
     private final Parameter meanParameter;
     private final Parameter precisionParameter;
     private final RandomField.WeightProvider weightProvider;
@@ -117,7 +117,7 @@ public double[] getMean() {
         return mean;
     }
 
-    SymmetricTriDiagonalMatrix getQ() {
+    protected SymmetricTriDiagonalMatrix getQ() {
         if (!qKnown) {
             double precision = precisionParameter.getParameterValue(0);
             Q.diagonal[0] = precision;