diff --git a/Figaro/META-INF/MANIFEST.MF b/Figaro/META-INF/MANIFEST.MF
index 2b33889f..52ab34a2 100644
--- a/Figaro/META-INF/MANIFEST.MF
+++ b/Figaro/META-INF/MANIFEST.MF
@@ -2,7 +2,7 @@ Manifest-Version: 1.0
 Bundle-ManifestVersion: 2
 Bundle-Name: Figaro
 Bundle-SymbolicName: com.cra.figaro
-Bundle-Version: 3.1.0
+Bundle-Version: 3.2.0
 Export-Package: com.cra.figaro.algorithm,
  com.cra.figaro.algorithm.decision,
  com.cra.figaro.algorithm.decision.index,
diff --git a/Figaro/figaro_build.properties b/Figaro/figaro_build.properties
index 9b3f6b46..cdb3799a 100644
--- a/Figaro/figaro_build.properties
+++ b/Figaro/figaro_build.properties
@@ -1 +1 @@
-version=3.1.0.0
+version=3.2.0.0
diff --git a/Figaro/src/main/scala/com/cra/figaro/algorithm/decision/DecisionVariableElimination.scala b/Figaro/src/main/scala/com/cra/figaro/algorithm/decision/DecisionVariableElimination.scala
index e9ce5c16..66c9f744 100644
--- a/Figaro/src/main/scala/com/cra/figaro/algorithm/decision/DecisionVariableElimination.scala
+++ b/Figaro/src/main/scala/com/cra/figaro/algorithm/decision/DecisionVariableElimination.scala
@@ -1,13 +1,13 @@
 /*
  * DecisionVariableElimination.scala
  * Variable elimination for Decisions algorithm.
- * 
+ *
  * Created By:      Brian Ruttenberg (bruttenberg@cra.com)
  * Creation Date:   Oct 1, 2012
- * 
+ *
  * Copyright 2013 Avrom J. Pfeffer and Charles River Analytics, Inc.
  * See http://www.cra.com or email figaro@cra.com for information.
- * 
+ *
  * See http://www.github.com/p2t2/figaro for a copy of the software license.
  */
 
@@ -28,7 +28,7 @@ import scala.language.existentials
 
 /* Trait only extends for double utilities. User needs to provide another trait or convert utilities to double
  * in order to use
- * 
+ *
  */
 /**
  * Trait for Decision based Variable Elimination. This implementation is hardcoded to use.
@@ -39,12 +39,12 @@ trait ProbabilisticVariableEliminationDecision extends VariableElimination[(Doub
    */
   /* Implementations must define this */
   def getUtilityNodes: List[Element[_]]
-  
+
   /**
    * Semiring for Decisions uses a sum-product-utility semiring.
    */
   override val semiring = SumProductUtilitySemiring()
-  
+
   /**
    * Makes a utility factor an element designated as a utility. This is factor of a tuple (Double, Double)
    * where the first value is 1.0 and the second is a possible utility of the element.
@@ -60,14 +60,14 @@ trait ProbabilisticVariableEliminationDecision extends VariableElimination[(Doub
   override def starterElements = getUtilityNodes ::: targetElements
 
   /**
-   * Create the factors for decision factors. Each factor is hardcoded as a tuple of (Double, Double), 
-   * where the first value is the probability and the second is the utility. 
-   */ 
+   * Create the factors for decision factors. Each factor is hardcoded as a tuple of (Double, Double),
+   * where the first value is the probability and the second is the utility.
+   */
   def getFactors(neededElements: List[Element[_]], targetElements: List[Element[_]], upper: Boolean = false): List[Factor[(Double, Double)]] = {
     if (debug) {
       println("Elements (other than utilities) appearing in factors and their ranges:")
-      for { element <- neededElements } { 
-        println(Variable(element).id + "(" + element.name.string + "@" + element.hashCode + ")" + ": " + element + ": " + Variable(element).range.mkString(",")) 
+      for { element <- neededElements } {
+        println(Variable(element).id + "(" + element.name.string + "@" + element.hashCode + ")" + ": " + element + ": " + Variable(element).range.mkString(","))
       }
     }
 
@@ -99,7 +99,7 @@ trait ProbabilisticVariableEliminationDecision extends VariableElimination[(Doub
       f.mapTo[(Double, Double)]((d: Double) => (d, 0.0), semiring.asInstanceOf[Semiring[(Double, Double)]])
     } else {
       if (f.variables.length > 1) throw new IllegalUtilityNodeException
-      
+
       val newF = f.mapTo[(Double, Double)]((d: Double) => (d, 0.0), semiring.asInstanceOf[Semiring[(Double, Double)]])
       for {i <- 0 until f.variables(0).range.size} {
         newF.set(List(i), (newF.get(List(i))._1, f.variables(0).range(i).asInstanceOf[Double]))
@@ -125,17 +125,17 @@ class ProbQueryVariableEliminationDecision[T, U](override val universe: Universe
   lazy val queryTargets = List(target)
 
   /**
-   *  The variable elimination eliminates all variables except on all decision nodes and their parents. 
+   *  The variable elimination eliminates all variables except on all decision nodes and their parents.
    *  Thus the target elements is both the decision element and the parent element.
    */
-  val targetElements = List(target, target.args(0)) 
+  val targetElements = List(target, target.args(0))
 
   def getUtilityNodes = utilityNodes
 
   private var finalFactors: Factor[(Double, Double)] = Factory.defaultFactor[(Double, Double)](List(), List(), semiring)
 
   /* Marginalizes the final factor using the semiring for decisions
-   * 
+   *
    */
   private def marginalizeToTarget(factor: Factor[(Double, Double)], target: Element[_]): Unit = {
     val unnormalizedTargetFactor = factor.marginalizeTo(semiring, Variable(target))
@@ -148,13 +148,13 @@ class ProbQueryVariableEliminationDecision[T, U](override val universe: Universe
   private def marginalize(resultFactor: Factor[(Double, Double)]) =
     queryTargets foreach (marginalizeToTarget(resultFactor, _))
 
-  private def makeResultFactor(factorsAfterElimination: Set[Factor[(Double, Double)]]): Factor[(Double, Double)] = {
+  private def makeResultFactor(factorsAfterElimination: MultiSet[Factor[(Double, Double)]]): Factor[(Double, Double)] = {
     // It is possible that there are no factors (this will happen if there are no decisions or utilities).
     // Therefore, we start with the unit factor and use foldLeft, instead of simply reducing the factorsAfterElimination.
     factorsAfterElimination.foldLeft(Factory.unit(semiring))(_.product(_))
   }
 
-  def finish(factorsAfterElimination: Set[Factor[(Double, Double)]], eliminationOrder: List[Variable[_]]) =
+  def finish(factorsAfterElimination: MultiSet[Factor[(Double, Double)]], eliminationOrder: List[Variable[_]]) =
     finalFactors = makeResultFactor(factorsAfterElimination)
 
   /**
@@ -176,7 +176,7 @@ class ProbQueryVariableEliminationDecision[T, U](override val universe: Universe
     (0.0 /: computeDistribution(target))(_ + get(_))
   }
 
-  /** 
+  /**
    *  Returns the computed utility of all parent/decision tuple values. For VE, these are not samples
    *  but the actual computed expected utility for all combinations of the parent and decision.
    */
@@ -194,14 +194,14 @@ class ProbQueryVariableEliminationDecision[T, U](override val universe: Universe
     // find the variables of the parents.
     val parentVariable = factor.variables.filterNot(_ == decisionVariable)(0)
 
-    // index of the decision variable     
+    // index of the decision variable
 
     val indexOfDecision = indices(factor.variables, decisionVariable)
     val indexOfParent = indices(factor.variables, parentVariable)
 
     for { indices <- factor.getIndices} {
 
-      /* for each index in the list of indices, strip out the decision variable index, 
+      /* for each index in the list of indices, strip out the decision variable index,
        * and retrieve the map entry for the parents. If the factor value is greater than
        * what is currently stored in the strategy map, replace the decision with the new one from the factor
        */
@@ -257,7 +257,7 @@ object DecisionVariableElimination {
   }
   /**
    * Create a decision variable elimination algorithm with  the given decision variables and indicated utility
-   * nodes and using the given dependent universes in the current default universe. Use the given dependent 
+   * nodes and using the given dependent universes in the current default universe. Use the given dependent
    * algorithm function to determine the algorithm to use to compute probability of evidence in each dependent universe.
    */
   def apply[T, U](
diff --git a/Figaro/src/main/scala/com/cra/figaro/algorithm/factored/FactoredAlgorithm.scala b/Figaro/src/main/scala/com/cra/figaro/algorithm/factored/FactoredAlgorithm.scala
index 0e2f327a..5747db53 100644
--- a/Figaro/src/main/scala/com/cra/figaro/algorithm/factored/FactoredAlgorithm.scala
+++ b/Figaro/src/main/scala/com/cra/figaro/algorithm/factored/FactoredAlgorithm.scala
@@ -73,7 +73,8 @@ trait FactoredAlgorithm[T] extends Algorithm {
         if (depth >= 0) {
           val includeThisElement = depth > LazyValues(element.universe).expandedDepth(element).getOrElse(-1)
           // Keeping track of what's been chased so far avoids infinite recursion
-          val toChase = element.universe.directlyUsedBy(element).toSet ++ dependentUniverseCoparents.getOrElse(element, Set()) -- chasedSoFar
+          val related = element.universe.directlyUsedBy(element).toSet ++ dependentUniverseCoparents.getOrElse(element, Set())
+          val toChase = related.filter(!chasedSoFar.contains(_))
           val rest = toChase.flatMap((elem: Element[_]) => chaseDown(elem, depth - 1, chasedSoFar ++ toChase))
           if (includeThisElement) rest + ((element, depth)); else rest
         } else Set()
diff --git a/Figaro/src/main/scala/com/cra/figaro/algorithm/factored/MPEVariableElimination.scala b/Figaro/src/main/scala/com/cra/figaro/algorithm/factored/MPEVariableElimination.scala
index a8e1f2c3..709c3194 100644
--- a/Figaro/src/main/scala/com/cra/figaro/algorithm/factored/MPEVariableElimination.scala
+++ b/Figaro/src/main/scala/com/cra/figaro/algorithm/factored/MPEVariableElimination.scala
@@ -1,13 +1,13 @@
 /*
  * MPEVariableElimination.scala
  * Variable elimination algorithm for computing most probable explanation.
- * 
+ *
  * Created By:      Avi Pfeffer (apfeffer@cra.com)
  * Creation Date:   Jan 1, 2009
- * 
+ *
  * Copyright 2013 Avrom J. Pfeffer and Charles River Analytics, Inc.
  * See http://www.cra.com or email figaro@cra.com for information.
- * 
+ *
  * See http://www.github.com/p2t2/figaro for a copy of the software license.
  */
 
@@ -20,10 +20,11 @@ import com.cra.figaro.algorithm.factored.factors._
 import com.cra.figaro.util
 import scala.collection.mutable.{ Set, Map }
 import com.cra.figaro.algorithm.lazyfactored._
+import com.cra.figaro.util.MultiSet
 
 /**
  * Variable elimination algorithm to compute the most probable explanation.
- * 
+ *
  * @param showTiming Produce timing information on steps of the algorithm
  */
 class MPEVariableElimination(override val universe: Universe)(
@@ -33,12 +34,12 @@ class MPEVariableElimination(override val universe: Universe)(
 
   override val comparator = Some((x: Double, y: Double) => x < y)
   override val semiring = MaxProductSemiring()
-  
+
   /*
    * We are trying to find a configuration of all the elements, so we must make them all starter elements for expansion.
    */
   override val starterElements = universe.activeElements
-  
+
   /**
    * Empty for MPE Algorithms.
    */
@@ -47,15 +48,15 @@ class MPEVariableElimination(override val universe: Universe)(
   private val maximizers: Map[Variable[_], Any] = Map()
 
   private def getMaximizer[T](variable: Variable[T]): T = maximizers(variable).asInstanceOf[variable.Value]
-  
+
   /*
    * Convert factors to use MaxProduct
    */
   override def getFactors(allElements: List[Element[_]], targetElements: List[Element[_]], upper: Boolean = false): List[Factor[Double]] = {
-    val factors = super.getFactors(allElements, targetElements, upper) 
+    val factors = super.getFactors(allElements, targetElements, upper)
     factors.map (_.mapTo(x => x, semiring))
   }
-  
+
   def mostLikelyValue[T](target: Element[T]): T = getMaximizer(Variable(target))
 
   private def backtrackOne[T](factor: Factor[_], variable: Variable[T]): Unit = {
@@ -64,7 +65,7 @@ class MPEVariableElimination(override val universe: Universe)(
     maximizers += variable -> factor.get(indices)
   }
 
-  def finish(factorsAfterElimination: Set[Factor[Double]], eliminationOrder: List[Variable[_]]): Unit =
+  def finish(factorsAfterElimination: MultiSet[Factor[Double]], eliminationOrder: List[Variable[_]]): Unit =
     for { (variable, factor) <- eliminationOrder.reverse.zip(recordingFactors) } { backtrackOne(factor, variable) }
 }
 
diff --git a/Figaro/src/main/scala/com/cra/figaro/algorithm/factored/SufficientStatisticsVariableElimination.scala b/Figaro/src/main/scala/com/cra/figaro/algorithm/factored/SufficientStatisticsVariableElimination.scala
index 3cf5786d..f37aaa41 100644
--- a/Figaro/src/main/scala/com/cra/figaro/algorithm/factored/SufficientStatisticsVariableElimination.scala
+++ b/Figaro/src/main/scala/com/cra/figaro/algorithm/factored/SufficientStatisticsVariableElimination.scala
@@ -1,13 +1,13 @@
 /*
  * SufficientStatisticsVariableElimination.scala
  * Variable elimination algorithm for sufficient statistics factors
- * 
+ *
  * Created By:      Michael Howard (mhoward@cra.com)
  * Creation Date:   Jun 1, 2013
- * 
+ *
  * Copyright 2013 Avrom J. Pfeffer and Charles River Analytics, Inc.
  * See http://www.cra.com or email figaro@cra.com for information.
- * 
+ *
  * See http://www.github.com/p2t2/figaro for a copy of the software license.
  */
 
@@ -19,12 +19,13 @@ import com.cra.figaro.language._
 import com.cra.figaro.algorithm.factored.factors._
 import scala.collection._
 import scala.collection.mutable.{ Map, Set }
+import com.cra.figaro.util.MultiSet
 
 /**
- * Variable elimination for sufficient statistics factors. 
+ * Variable elimination for sufficient statistics factors.
  * The final factor resulting from variable elimination contains a mapping of parameters to sufficient statistics vectors
  * which can be used to maximize parameter values.
- * 
+ *
  * @param parameterMap A map of parameters to their sufficient statistics.
  */
 class SufficientStatisticsVariableElimination(
@@ -47,14 +48,14 @@ class SufficientStatisticsVariableElimination(
   }
 
   /**
-   *  Particular implementations of probability of evidence algorithms must define the following method. 
+   *  Particular implementations of probability of evidence algorithms must define the following method.
    */
   def getFactors(neededElements: List[Element[_]], targetElements: List[Element[_]], upper: Boolean = false): List[Factor[(Double, mutable.Map[Parameter[_], Seq[Double]])]] = {
     val allElements = neededElements.filter(p => p.isInstanceOf[Parameter[_]] == false)
     if (debug) {
       println("Elements appearing in factors and their ranges:")
-      for { element <- allElements } { 
-        println(Variable(element).id + "(" + element.name.string + "@" + element.hashCode + ")" + ": " + element + ": " + Variable(element).range.mkString(",")) 
+      for { element <- allElements } {
+        println(Variable(element).id + "(" + element.name.string + "@" + element.hashCode + ")" + ": " + element + ": " + Variable(element).range.mkString(","))
       }
     }
 
@@ -70,12 +71,12 @@ class SufficientStatisticsVariableElimination(
    * Empty for this algorithm.
    */
   val targetElements = List[Element[_]]()
-  
+
   override def starterElements = universe.conditionedElements ++ universe.constrainedElements
 
   private var result: (Double, Map[Parameter[_], Seq[Double]]) = _
 
-  def finish(factorsAfterElimination: Set[Factor[(Double, Map[Parameter[_], Seq[Double]])]], eliminationOrder: List[Variable[_]]): Unit = {
+  def finish(factorsAfterElimination: MultiSet[Factor[(Double, Map[Parameter[_], Seq[Double]])]], eliminationOrder: List[Variable[_]]): Unit = {
     // It is possible that there are no factors (this will happen if there is no evidence).
     // Therefore, we start with the unit factor and use foldLeft, instead of simply reducing the factorsAfterElimination.
     val finalFactor = factorsAfterElimination.foldLeft(Factory.unit(semiring))(_.product(_))
@@ -86,14 +87,14 @@ class SufficientStatisticsVariableElimination(
   }
 
   /**
-   * Returns a mapping of parameters to sufficient statistics resulting from 
+   * Returns a mapping of parameters to sufficient statistics resulting from
    * elimination of the factors.
    */
   def getSufficientStatisticsForAllParameters = { result._2.toMap }
 
   val semiring = SufficientStatisticsSemiring(parameterMap)
 
-  override def cleanUp() = { 
+  override def cleanUp() = {
     statFactor.removeFactors
     super.cleanUp()
   }
diff --git a/Figaro/src/main/scala/com/cra/figaro/algorithm/factored/VariableElimination.scala b/Figaro/src/main/scala/com/cra/figaro/algorithm/factored/VariableElimination.scala
index a1d49c97..ec030306 100644
--- a/Figaro/src/main/scala/com/cra/figaro/algorithm/factored/VariableElimination.scala
+++ b/Figaro/src/main/scala/com/cra/figaro/algorithm/factored/VariableElimination.scala
@@ -66,13 +66,17 @@ trait VariableElimination[T] extends FactoredAlgorithm[T] with OneTime {
 
   // The first element of FactorMap is the complete set of factors.
   // The second element maps variables to the factors mentioning that variable.
-  private type FactorMap[T] = Map[Variable[_], Set[Factor[T]]]
+  // The previous implementation used sets, but that resulted in bugs where an identical factor appeared more than once.
+  // The new implementation uses multisets.
+  private type FactorMap[T] = Map[Variable[_], MultiSet[Factor[T]]]
 
+  // Add a factor to the list, even if it appears already.
   private def addFactor[T](factor: Factor[T], map: FactorMap[T]): Unit =
-    factor.variables foreach (v => map += v -> (map.getOrElse(v, Set()) + factor))
+    factor.variables foreach (v => map += v -> (map.getOrElse(v, HashMultiSet()).addOne(factor)))
 
+  // Remove one instance of the factor from the list.
   private def removeFactor[T](factor: Factor[T], map: FactorMap[T]): Unit =
-    factor.variables foreach (v => map += v -> (map.getOrElse(v, Set()) - factor))
+    factor.variables foreach (v => map += v -> (map.getOrElse(v, HashMultiSet()).removeOne(factor)))
 
   protected def initialFactorMap(factors: Traversable[Factor[T]]): FactorMap[T] = {
     val map: FactorMap[T] = Map()
@@ -95,7 +99,7 @@ trait VariableElimination[T] extends FactoredAlgorithm[T] with OneTime {
 
   private def eliminate(
     variable: Variable[_],
-    factors: Set[Factor[T]],
+    factors: MultiSet[Factor[T]],
     map: FactorMap[T]): Unit = {
     val varFactors = map(variable)
     if (debug) {
@@ -106,23 +110,23 @@ trait VariableElimination[T] extends FactoredAlgorithm[T] with OneTime {
     if (varFactors nonEmpty) {
       val productFactor = varFactors reduceLeft (_.product(_))
       val resultFactor = productFactor.sumOver(variable)
-      varFactors foreach (removeFactor(_, map))
-      addFactor(resultFactor, map)
+      if (debug) println("Result factor\n" + resultFactor.toReadableString)
       comparator match {
         case None => ()
         case Some(recorder) => recordingFactors ::= productFactor.recordArgMax(variable, recorder)
       }
+      varFactors.foreach(factors.removeOne(_))
+      factors.addOne(resultFactor)
+      varFactors.foreach(removeFactor(_, map))
       map -= variable
-      factors --= varFactors
-      if (debug) println("Result factor\n" + resultFactor.toReadableString)
-      factors += resultFactor
+      addFactor(resultFactor, map)
     }
   }
 
   protected def eliminateInOrder(
     order: List[Variable[_]],
-    factors: Set[Factor[T]],
-    map: FactorMap[T]): Set[Factor[T]] =
+    factors: MultiSet[Factor[T]],
+    map: FactorMap[T]): MultiSet[Factor[T]] =
     order match {
       case Nil =>
         factors
@@ -147,7 +151,7 @@ trait VariableElimination[T] extends FactoredAlgorithm[T] with OneTime {
     }
     val (_, order) = optionallyShowTiming(VariableElimination.eliminationOrder(allFactors, targetVariables), "Computing elimination order")
     val factorsAfterElimination =
-      optionallyShowTiming(eliminateInOrder(order, Set(allFactors: _*), initialFactorMap(allFactors)), "Elimination")
+      optionallyShowTiming(eliminateInOrder(order, HashMultiSet(allFactors: _*), initialFactorMap(allFactors)), "Elimination")
     if (debug) println("*****************")
     if (debug) factorsAfterElimination foreach (f => println(f.toReadableString))
     optionallyShowTiming(finish(factorsAfterElimination, order), "Finalizing")
@@ -159,7 +163,7 @@ trait VariableElimination[T] extends FactoredAlgorithm[T] with OneTime {
   /**
    * All implementation of variable elimination must specify what to do after variables have been eliminated.
    */
-  def finish(factorsAfterElimination: Set[Factor[T]], eliminationOrder: List[Variable[_]]): Unit
+  def finish(factorsAfterElimination: MultiSet[Factor[T]], eliminationOrder: List[Variable[_]]): Unit
 
   def run() = ve()
 
@@ -211,13 +215,13 @@ class ProbQueryVariableElimination(override val universe: Universe, targets: Ele
   private def marginalize(resultFactor: Factor[Double]) =
     targets foreach (marginalizeToTarget(resultFactor, _))
 
-  private def makeResultFactor(factorsAfterElimination: Set[Factor[Double]]): Factor[Double] = {
+  private def makeResultFactor(factorsAfterElimination: MultiSet[Factor[Double]]): Factor[Double] = {
     // It is possible that there are no factors (this will happen if there are  no queries or evidence).
     // Therefore, we start with the unit factor and use foldLeft, instead of simply reducing the factorsAfterElimination.
     factorsAfterElimination.foldLeft(Factory.unit(semiring))(_.product(_))
   }
 
-  def finish(factorsAfterElimination: Set[Factor[Double]], eliminationOrder: List[Variable[_]]) =
+  def finish(factorsAfterElimination: MultiSet[Factor[Double]], eliminationOrder: List[Variable[_]]) =
     marginalize(makeResultFactor(factorsAfterElimination))
 
   /**
diff --git a/Figaro/src/main/scala/com/cra/figaro/algorithm/factored/factors/factory/ComplexFactory.scala b/Figaro/src/main/scala/com/cra/figaro/algorithm/factored/factors/factory/ComplexFactory.scala
index 679b8b35..fb189b46 100644
--- a/Figaro/src/main/scala/com/cra/figaro/algorithm/factored/factors/factory/ComplexFactory.scala
+++ b/Figaro/src/main/scala/com/cra/figaro/algorithm/factored/factors/factory/ComplexFactory.scala
@@ -1,13 +1,13 @@
 /*
  * ComplexFactory.scala
  * Description needed
- * 
+ *
  * Created By:      Glenn Takata (gtakata@cra.com)
  * Creation Date:   Dec 15, 2014
- * 
+ *
  * Copyright 2014 Avrom J. Pfeffer and Charles River Analytics, Inc.
  * See http://www.cra.com or email figaro@cra.com for information.
- * 
+ *
  * See http://www.github.com/p2t2/figaro for a copy of the software license.
  */
 
@@ -53,7 +53,7 @@ object ComplexFactory {
         selectedFactors.flatten
     }
   }
-  
+
   /**
    * Factor constructor for a MultiValuedReferenceElement
    */
@@ -171,7 +171,7 @@ object ComplexFactory {
     }
     List(factor)
   }
-  
+
   /**
    * Factor constructor for a FoldLeft Element
    */
diff --git a/Figaro/src/main/scala/com/cra/figaro/algorithm/learning/GeneralizedEM.scala b/Figaro/src/main/scala/com/cra/figaro/algorithm/learning/GeneralizedEM.scala
index 27d0df86..0d9c6e21 100644
--- a/Figaro/src/main/scala/com/cra/figaro/algorithm/learning/GeneralizedEM.scala
+++ b/Figaro/src/main/scala/com/cra/figaro/algorithm/learning/GeneralizedEM.scala
@@ -83,21 +83,23 @@ trait ExpectationMaximization extends Algorithm with ParameterLearner {
  */
 trait OnlineExpectationMaximization extends Online with ExpectationMaximization {
 
+  override def doStart = {}
+
   protected var lastIterationStatistics: Map[Parameter[_], Seq[Double]] = Map(targetParameters.map(p => p -> p.zeroSufficientStatistics): _*)
   override val initial: Universe
   override val transition: Function0[Universe]
   protected var currentUniverse: Universe = initial
 
-  private def updateStatistics(newStatistics:Map[Parameter[_], Seq[Double]]): Map[Parameter[_], Seq[Double]] = {
+  private def updateStatistics(newStatistics: Map[Parameter[_], Seq[Double]]): Map[Parameter[_], Seq[Double]] = {
     Map((for (p <- paramMap.keys) yield {
       val updatedStatistics = (lastIterationStatistics(p) zip newStatistics(p)).map((pair: (Double, Double)) => pair._1 + pair._2)
       (p, updatedStatistics)
-    }).toSeq:_*)
-  } 
-  
- /**
- * Observe new evidence and perform one expectation step and one maximization step
- */
+    }).toSeq: _*)
+  }
+
+  /**
+   * Observe new evidence and perform one expectation step and one maximization step
+   */
   def update(evidence: Seq[NamedEvidence[_]] = Seq()): Unit = {
     currentUniverse = transition()
     currentUniverse.assertEvidence(evidence)
@@ -124,13 +126,12 @@ class ExpectationMaximizationWithFactors(val universe: Universe, val targetParam
 
 }
 
-
 /**
  * An online EM algorithm which learns parameters using a factored algorithm
  */
 class OnlineExpectationMaximizationWithFactors(override val initial: Universe, override val transition: Function0[Universe], val targetParameters: Parameter[_]*)(val terminationCriteria: () => EMTerminationCriteria)
   extends OnlineExpectationMaximization {
-  
+
   def doExpectationStep = {
     val algorithm = SufficientStatisticsVariableElimination(paramMap)(currentUniverse)
     algorithm.start
@@ -141,12 +142,12 @@ class OnlineExpectationMaximizationWithFactors(override val initial: Universe, o
   }
 }
 
-
 /**
  * An EM algorithm which learns parameters using an inference algorithm provided as an argument
  */
 class GeneralizedEM(inferenceAlgorithmConstructor: Seq[Element[_]] => Universe => ProbQueryAlgorithm with OneTime, val universe: Universe, val targetParameters: Parameter[_]*)(val terminationCriteria: () => EMTerminationCriteria) extends ExpectationMaximization {
 
+  //Dependent universe doesn't work the same way.
   protected def doExpectationStep(): Map[Parameter[_], Seq[Double]] = {
     val inferenceTargets =
       universe.activeElements.filter(_.isInstanceOf[Parameterized[_]]).map(_.asInstanceOf[Parameterized[_]])
@@ -161,6 +162,7 @@ class GeneralizedEM(inferenceAlgorithmConstructor: Seq[Element[_]] => Universe =
       for {
         target <- universe.directlyUsedBy(parameter)
       } {
+
         val t: Parameterized[target.Value] = target.asInstanceOf[Parameterized[target.Value]]
         val distribution: Stream[(Double, target.Value)] = algorithm.distribution(t)
         val newStats = t.distributionToStatistics(parameter, distribution)
@@ -179,21 +181,27 @@ class GeneralizedEM(inferenceAlgorithmConstructor: Seq[Element[_]] => Universe =
  */
 class GeneralizedOnlineEM(inferenceAlgorithmConstructor: Seq[Element[_]] => Universe => ProbQueryAlgorithm with OneTime, override val initial: Universe, override val transition: Function0[Universe], val targetParameters: Parameter[_]*)(val terminationCriteria: () => EMTerminationCriteria) extends OnlineExpectationMaximization {
 
-  
+  protected def usesParameter(l: List[Element[_]]): Map[Parameter[_], Iterable[Parameterized[_]]] = {
+    (l.map { x => x match { case p: Parameterized[_] => { p -> p.parameters.head } } }).groupBy(_._2).mapValues(_.map(_._1))
+  }
+
   protected def doExpectationStep(): Map[Parameter[_], Seq[Double]] = {
     val inferenceTargets =
       currentUniverse.activeElements.filter(_.isInstanceOf[Parameterized[_]]).map(_.asInstanceOf[Parameterized[_]])
 
     val algorithm = inferenceAlgorithmConstructor(inferenceTargets)(currentUniverse)
     algorithm.start()
-
+    //println("universe: " + currentUniverse.hashCode)
     var result: Map[Parameter[_], Seq[Double]] = Map()
 
+    val uses = usesParameter(inferenceTargets)
+    println("built map")
     for { parameter <- targetParameters } {
       var stats = parameter.zeroSufficientStatistics
       for {
-        target <- currentUniverse.directlyUsedBy(parameter)
+        target <- uses(parameter)
       } {
+        println("found used by...")
         val t: Parameterized[target.Value] = target.asInstanceOf[Parameterized[target.Value]]
         val distribution: Stream[(Double, target.Value)] = algorithm.distribution(t)
         val newStats = t.distributionToStatistics(parameter, distribution)
@@ -208,7 +216,7 @@ class GeneralizedOnlineEM(inferenceAlgorithmConstructor: Seq[Element[_]] => Univ
 }
 
 object EMWithBP {
-  
+
   private val defaultBPIterations = 10
 
   def online(transition: () => Universe, p: Parameter[_]*)(implicit universe: Universe) = {
@@ -270,9 +278,9 @@ object EMWithBP {
 }
 
 object EMWithImportance {
-  
+
   private val defaultImportanceParticles = 100000
-  
+
   private def makeImportance(numParticles: Int, targets: Seq[Element[_]])(universe: Universe) = {
     Importance(numParticles, targets: _*)(universe)
   }
@@ -342,7 +350,7 @@ object EMWithImportance {
 object EMWithMH {
 
   private val defaultMHParticles = 100000
-  
+
   private def makeImportance(numParticles: Int, targets: Seq[Element[_]])(universe: Universe) = {
     Importance(numParticles, targets: _*)(universe)
   }
diff --git a/Figaro/src/main/scala/com/cra/figaro/algorithm/sampling/Forward.scala b/Figaro/src/main/scala/com/cra/figaro/algorithm/sampling/Forward.scala
index 790cff64..e1b73a16 100644
--- a/Figaro/src/main/scala/com/cra/figaro/algorithm/sampling/Forward.scala
+++ b/Figaro/src/main/scala/com/cra/figaro/algorithm/sampling/Forward.scala
@@ -80,7 +80,7 @@ object Forward {
             var argsRemaining = initialArgs 
             while (!argsRemaining.isEmpty) {
               state1 = sampleInState(argsRemaining.head, state1, universe, useObservation)
-              val newArgs = element.args.toSet -- initialArgs
+              val newArgs = element.args.filter(!initialArgs.contains(_))
               initialArgs = initialArgs ++ newArgs
               argsRemaining = argsRemaining.tail ++ newArgs
             }
diff --git a/Figaro/src/main/scala/com/cra/figaro/algorithm/sampling/Importance.scala b/Figaro/src/main/scala/com/cra/figaro/algorithm/sampling/Importance.scala
index d5760dd1..b70ef303 100644
--- a/Figaro/src/main/scala/com/cra/figaro/algorithm/sampling/Importance.scala
+++ b/Figaro/src/main/scala/com/cra/figaro/algorithm/sampling/Importance.scala
@@ -251,8 +251,9 @@ abstract class Importance(universe: Universe, targets: Element[_]*)
   @tailrec
   private def sampleArgs(element: Element[_], state: State, args: Set[Element[_]]): Unit = {
     args foreach (sampleOne(state, _, None))
-    val newArgs = Set[Element[_]](element.args:_*) -- state.assigned
-    if (newArgs.nonEmpty) sampleArgs(element, state, newArgs) else return
+    val newArgs = element.args.filter(!state.assigned.contains(_))
+    if (newArgs.nonEmpty)
+      sampleArgs(element, state, Set[Element[_]](newArgs: _*))
   }
 
    /**
diff --git a/Figaro/src/main/scala/com/cra/figaro/algorithm/sampling/MetropolisHastings.scala b/Figaro/src/main/scala/com/cra/figaro/algorithm/sampling/MetropolisHastings.scala
index 0f5189f2..e727243a 100644
--- a/Figaro/src/main/scala/com/cra/figaro/algorithm/sampling/MetropolisHastings.scala
+++ b/Figaro/src/main/scala/com/cra/figaro/algorithm/sampling/MetropolisHastings.scala
@@ -217,37 +217,26 @@ abstract class MetropolisHastings(universe: Universe, proposalScheme: ProposalSc
    * is not empty, we recursively update the intersecting elements. Once those updates are completed,
    * we update an element and move on to the next element in the set.
    */
-  protected def updateMany[T](state: State, toUpdate: Set[Element[_]]): State = {
-    var returnState = state
-    var updatesLeft = toUpdate
-    while (!updatesLeft.isEmpty) {
-      // Check the intersection of an element's arguments with the updates that still need to occur
-      var argsRemaining = universe.uses(updatesLeft.head).intersect(updatesLeft)
-      while (!argsRemaining.isEmpty) {
-        // update the element's arguments first
-        returnState = updateManyHelper(returnState, argsRemaining.toSet)
-        argsRemaining = argsRemaining.tail
+  @tailrec
+  private def updateMany(state: State, currentStack: List[Element[_]], currentArgs: Set[Element[_]], updateQ: Set[Element[_]]): State = {
+
+    if (currentStack.isEmpty && currentArgs.isEmpty && updateQ.isEmpty) state
+
+    else if (currentStack.isEmpty && currentArgs.isEmpty && updateQ.nonEmpty) { 
+      val argsRemaining = universe.uses(updateQ.head).intersect(updateQ.tail)
+      updateMany(state, List(updateQ.head), argsRemaining.toSet, updateQ.tail -- argsRemaining)
+    } else if (currentStack.nonEmpty && currentArgs.isEmpty) { 
+      val newState = updateOne(state, currentStack.head)
+      updateMany(newState, currentStack.tail, currentArgs, updateQ)
+    } else {
+      val argsRemaining = universe.uses(currentArgs.head).intersect(currentArgs.tail)
+      if (argsRemaining.isEmpty) {
+        val newState = updateOne(state, currentArgs.head)
+        updateMany(newState, currentStack, currentArgs.tail, updateQ)
+      } else {
+        updateMany(state, currentArgs.head +: currentStack, currentArgs.tail, updateQ)
       }
-      // once the args are updated, update this element
-      returnState = updateOne(returnState, updatesLeft.head)
-      updatesLeft = updatesLeft.tail
     }
-    returnState
-  }
-
-  /*
-   * A recursive function to work in conjunction with updateMany to check the order of the element
-   * updates.
-   */
-  @tailrec
-  private def updateManyHelper(state: State, toUpdate: Set[Element[_]]): State = {
-    var returnState = state
-    var updatesLeft = toUpdate
-    var argsRemaining = universe.uses(updatesLeft.head).intersect(updatesLeft)
-    if (argsRemaining.isEmpty) {
-      returnState = updateOne(returnState, updatesLeft.head)
-      returnState
-    } else { updateManyHelper(returnState, argsRemaining.toSet) }
   }
 
   /*
@@ -257,7 +246,7 @@ abstract class MetropolisHastings(universe: Universe, proposalScheme: ProposalSc
   protected def proposeAndUpdate(): State = {
     val state1 = runScheme()
     val updatesNeeded = state1.oldValues.keySet flatMap (elem => universe.usedBy(elem))
-    updateMany(state1, updatesNeeded.toSet)
+    updateMany(state1, List(), Set(), updatesNeeded.toSet)
   }
 
   protected var dissatisfied: Set[Element[_]] = _
@@ -365,7 +354,7 @@ abstract class MetropolisHastings(universe: Universe, proposalScheme: ProposalSc
 
   protected def doInitialize(): Unit = {
     // Need to prime the universe to make sure all elements have a generated value
-    Forward(true)(universe)
+    Forward(false)(universe)
     initConstrainedValues()
     dissatisfied = universe.conditionedElements.toSet filter (!_.conditionSatisfied)
     for { i <- 1 to burnIn } mhStep()
diff --git a/Figaro/src/main/scala/com/cra/figaro/algorithm/sampling/MetropolisHastingsAnnealer.scala b/Figaro/src/main/scala/com/cra/figaro/algorithm/sampling/MetropolisHastingsAnnealer.scala
index d347f256..ed361e91 100644
--- a/Figaro/src/main/scala/com/cra/figaro/algorithm/sampling/MetropolisHastingsAnnealer.scala
+++ b/Figaro/src/main/scala/com/cra/figaro/algorithm/sampling/MetropolisHastingsAnnealer.scala
@@ -104,8 +104,7 @@ abstract class MetropolisHastingsAnnealer(universe: Universe, proposalScheme: Pr
 
     if (dissatisfied.isEmpty) {
       sampleCount += 1
-      val toUpdate = if (currentEnergy > bestEnergy) {
-        allLastUpdates.clear
+      val toUpdate = if (currentEnergy >= bestEnergy) {
         saveState
       } else Map[Element[_], Any]()
       (true, toUpdate)
@@ -115,7 +114,7 @@ abstract class MetropolisHastingsAnnealer(universe: Universe, proposalScheme: Pr
   }
 
   override def doInitialize(): Unit = {
-    Forward(true)(universe)
+    Forward(false)(universe)
     initConstrainedValues()
     dissatisfied = universe.conditionedElements.toSet filter (!_.conditionSatisfied)
     currentEnergy = universe.constrainedElements.map(_.constraintValue).sum
@@ -123,6 +122,8 @@ abstract class MetropolisHastingsAnnealer(universe: Universe, proposalScheme: Pr
       val nextState = mhStep()
       currentEnergy += nextState.modelProb
     }
+    initUpdates()    
+    if (dissatisfied.nonEmpty) bestEnergy = Double.MinValue else bestEnergy = currentEnergy 
   }
 
   def mostLikelyValue[T](target: Element[T]): T = {
diff --git a/Figaro/src/main/scala/com/cra/figaro/experimental/structured/Problem.scala b/Figaro/src/main/scala/com/cra/figaro/experimental/structured/Problem.scala
index de6d2278..0670c9d2 100644
--- a/Figaro/src/main/scala/com/cra/figaro/experimental/structured/Problem.scala
+++ b/Figaro/src/main/scala/com/cra/figaro/experimental/structured/Problem.scala
@@ -60,6 +60,15 @@ class Problem(val collection: ComponentCollection, targets: List[Element[_]] = L
       contains(component.problem) & !targets.contains(component.element)
     }
   }
+  
+  /**
+   * Determines if a variable is in scope outside of this problem
+   */
+  def global(variable: Variable[_]): Boolean = {
+    !collection.intermediates.contains(variable) &&
+    !contains(collection.variableToComponent(variable).problem)
+  }
+  
 
   /**
    * Determines if this problem contains the given problem.
diff --git a/Figaro/src/main/scala/com/cra/figaro/experimental/structured/ProblemComponent.scala b/Figaro/src/main/scala/com/cra/figaro/experimental/structured/ProblemComponent.scala
index 91467dab..3fbb0158 100644
--- a/Figaro/src/main/scala/com/cra/figaro/experimental/structured/ProblemComponent.scala
+++ b/Figaro/src/main/scala/com/cra/figaro/experimental/structured/ProblemComponent.scala
@@ -189,7 +189,7 @@ extends ExpandableComponent[ParentValue, Value](problem, chain.parent, chain) {
         (parentValue, subproblem) <- subproblems
         factor <- subproblem.solution
       } yield Factory.replaceVariable(factor, problem.collection(subproblem.target).variable, actualSubproblemVariables(parentValue))
-    nonConstraintFactors = subproblems.values.toList.flatMap(_.solution) ::: nonConstraintFactors
+    nonConstraintFactors = subproblemFactors.toList ::: nonConstraintFactors
   }
 
   /*
diff --git a/Figaro/src/main/scala/com/cra/figaro/experimental/structured/factory/ChainFactory.scala b/Figaro/src/main/scala/com/cra/figaro/experimental/structured/factory/ChainFactory.scala
index 8750c30f..9e8b1d5d 100644
--- a/Figaro/src/main/scala/com/cra/figaro/experimental/structured/factory/ChainFactory.scala
+++ b/Figaro/src/main/scala/com/cra/figaro/experimental/structured/factory/ChainFactory.scala
@@ -43,7 +43,7 @@ object ChainFactory {
         val nestedProblem = cc.expansions((chain.chainFunction, parentVal.value))
         val outcomeElem = nestedProblem.target.asInstanceOf[Element[U]]
         val formalVar = Factory.getVariable(cc, outcomeElem)
-        val actualVar = if (nestedProblem.internal(formalVar)) Factory.makeVariable(cc, formalVar.valueSet) else formalVar
+        val actualVar = if (!nestedProblem.global(formalVar)) Factory.makeVariable(cc, formalVar.valueSet) else formalVar
         chainComp.actualSubproblemVariables += parentVal.value -> actualVar
         List(Factory.makeConditionalSelector(pairVar, parentVal, actualVar))
       }
diff --git a/Figaro/src/main/scala/com/cra/figaro/experimental/structured/solver/VEBPChooser.scala b/Figaro/src/main/scala/com/cra/figaro/experimental/structured/solver/VEBPChooser.scala
index 5acd9d47..77af890c 100644
--- a/Figaro/src/main/scala/com/cra/figaro/experimental/structured/solver/VEBPChooser.scala
+++ b/Figaro/src/main/scala/com/cra/figaro/experimental/structured/solver/VEBPChooser.scala
@@ -15,6 +15,8 @@ package com.cra.figaro.experimental.structured.solver
 import com.cra.figaro.algorithm.factored.VariableElimination
 import com.cra.figaro.algorithm.factored.factors._
 import com.cra.figaro.experimental.structured.Problem
+import com.cra.figaro.util.HashMultiSet
+import com.cra.figaro.util.MultiSet
 
 /*
  * The score threshold is the threshold to determine whether VE or BP will be used.
@@ -33,7 +35,7 @@ extends VariableElimination[Double] {
     } else {
 //      println("Choosing VE")
       // Since we've already computed the order, we don't call doElimination but only do the steps after computing the order
-      val factorsAfterElimination = eliminateInOrder(order, scala.collection.mutable.Set(factors: _*), initialFactorMap(factors))
+      val factorsAfterElimination = eliminateInOrder(order, HashMultiSet(factors: _*), initialFactorMap(factors))
       finish(factorsAfterElimination, order)
     }
     result
@@ -43,7 +45,7 @@ extends VariableElimination[Double] {
 
   private var result: List[Factor[Double]] = _
 
-  def finish(factorsAfterElimination: scala.collection.mutable.Set[Factor[Double]], eliminationOrder: List[Variable[_]]): Unit = {
+  def finish(factorsAfterElimination: MultiSet[Factor[Double]], eliminationOrder: List[Variable[_]]): Unit = {
     result = factorsAfterElimination.toList
   }
 
diff --git a/Figaro/src/main/scala/com/cra/figaro/experimental/structured/solver/VESolver.scala b/Figaro/src/main/scala/com/cra/figaro/experimental/structured/solver/VESolver.scala
index b84a7044..27a7a006 100644
--- a/Figaro/src/main/scala/com/cra/figaro/experimental/structured/solver/VESolver.scala
+++ b/Figaro/src/main/scala/com/cra/figaro/experimental/structured/solver/VESolver.scala
@@ -14,6 +14,7 @@ package com.cra.figaro.experimental.structured.solver
 
 import com.cra.figaro.algorithm.factored.factors._
 import com.cra.figaro.experimental.structured.Problem
+import com.cra.figaro.util.MultiSet
 
 private[figaro] class VESolver(problem: Problem, toEliminate: Set[Variable[_]], toPreserve: Set[Variable[_]], factors: List[Factor[Double]])
 extends com.cra.figaro.algorithm.factored.VariableElimination[Double] {
@@ -22,13 +23,13 @@ extends com.cra.figaro.algorithm.factored.VariableElimination[Double] {
     result
   }
 
-  debug = true
+  debug = false
 
   val semiring: Semiring[Double] = SumProductSemiring()
 
   private var result: List[Factor[Double]] = _
 
-  def finish(factorsAfterElimination: scala.collection.mutable.Set[Factor[Double]], eliminationOrder: List[Variable[_]]): Unit = {
+  def finish(factorsAfterElimination: MultiSet[Factor[Double]], eliminationOrder: List[Variable[_]]): Unit = {
     result = factorsAfterElimination.toList
   }
 
diff --git a/Figaro/src/test/scala/com/cra/figaro/test/algorithm/factored/VETest.scala b/Figaro/src/test/scala/com/cra/figaro/test/algorithm/factored/VETest.scala
index b2ed3474..f41ad157 100644
--- a/Figaro/src/test/scala/com/cra/figaro/test/algorithm/factored/VETest.scala
+++ b/Figaro/src/test/scala/com/cra/figaro/test/algorithm/factored/VETest.scala
@@ -1,13 +1,13 @@
 /*
- * VETest.scala  
+ * VETest.scala
  * Variable elimination tests.
- * 
+ *
  * Created By:      Avi Pfeffer (apfeffer@cra.com)
  * Creation Date:   Jan 1, 2009
- * 
+ *
  * Copyright 2013 Avrom J. Pfeffer and Charles River Analytics, Inc.
  * See http://www.cra.com or email figaro@cra.com for information.
- * 
+ *
  * See http://www.github.com/p2t2/figaro for a copy of the software license.
  */
 
@@ -306,7 +306,7 @@ class VETest extends WordSpec with Matchers {
       val a = If(f, Select(0.3 -> 1, 0.7 -> 2), Constant(2))
       test(f, (b: Boolean) => b, 0.6)
     }
-    
+
     "on a different universe from the current universe, produce the correct result" in {
       val u1 = Universe.createNew()
       val u = Select(0.25 -> 0.3, 0.25 -> 0.5, 0.25 -> 0.7, 0.25 -> 0.9)
@@ -381,7 +381,7 @@ class VETest extends WordSpec with Matchers {
       ve.probability(y, true) should be(((0.1 * 0.2 + 0.9 * 0.2) / (0.1 * 0.2 + 0.9 * 0.2 + 0.9 * 0.8)) +- 0.0000000001)
       ve.kill
     }
-    
+
   }
 
   "MPEVariableElimination" should {
@@ -407,7 +407,7 @@ class VETest extends WordSpec with Matchers {
       alg.kill
     }
   }
-  
+
   def test[T](target: Element[T], predicate: T => Boolean, prob: Double) {
     val tolerance = 0.0000001
     val algorithm = VariableElimination(target)
diff --git a/Figaro/src/test/scala/com/cra/figaro/test/algorithm/learning/OnlineEMTest.scala b/Figaro/src/test/scala/com/cra/figaro/test/algorithm/learning/OnlineEMTest.scala
index f1d48104..205e31fa 100644
--- a/Figaro/src/test/scala/com/cra/figaro/test/algorithm/learning/OnlineEMTest.scala
+++ b/Figaro/src/test/scala/com/cra/figaro/test/algorithm/learning/OnlineEMTest.scala
@@ -29,7 +29,239 @@ import scala.math.abs
 import java.io._
 import com.cra.figaro.test.tags.NonDeterministic
 
-class OnlineEMTest extends WordSpec with PrivateMethodTester with Matchers { 
+class OnlineEMTest extends WordSpec with PrivateMethodTester with Matchers {
+  
+    "Online expectation maximization with BP" when
+    {
+
+      "used to estimate a Beta parameter" should
+        {
+
+          "detect bias after a large enough number of trials" in
+            {
+              val initial = Universe.createNew
+              val b = Beta(2, 2)
+              def transition = () => {
+                val next = Universe.createNew
+                val f = Flip(b)("f",next)
+                next
+              }
+              val algorithm = EMWithBP.online(transition, b)(initial)
+              algorithm.start()
+              for (i <- 1 to 7) {
+                algorithm.update(List(NamedEvidence("f", Observation(true))))
+              }
+
+              for (i <- 1 to 3) {
+                algorithm.update(List(NamedEvidence("f", Observation(false))))
+              }
+
+
+              val result = b.MAPValue
+              algorithm.kill
+              result should be(0.6666 +- 0.01)
+
+            }
+
+          "take the prior concentration parameters into account" in
+            {
+              val initial = Universe.createNew
+              val b = Beta(3.0, 7.0)
+              def transition = () => {
+                val next = Universe.createNew
+                val f = Flip(b)("f",next)
+                next
+              }
+              val algorithm = EMWithBP.online(transition, b)(initial)
+              algorithm.start()
+              for (i <- 1 to 7) {
+                algorithm.update(List(NamedEvidence("f", Observation(true))))
+              }
+              for (i <- 1 to 3) {
+                algorithm.update(List(NamedEvidence("f", Observation(false))))
+              }
+              val result = b.MAPValue
+              algorithm.kill
+              result should be(0.50 +- 0.01)
+            }
+
+          "learn the bias from observations of binomial elements" in {
+              val initial = Universe.createNew
+              val b = Beta(2, 2)
+              def transition = () => {
+                val next = Universe.createNew
+                val f = Binomial(5,b)("binomial",next)
+                next
+              }
+              val algorithm = EMWithBP.online(transition, b)(initial)
+              algorithm.start()
+              algorithm.update(List(NamedEvidence("binomial", Observation(4))))
+              algorithm.update(List(NamedEvidence("binomial", Observation(3))))
+              val result = b.MAPValue
+              algorithm.kill
+              result should be(0.6666 +- 0.01)
+          }
+
+          "correctly use a uniform prior" in {
+              val initial = Universe.createNew
+              val b = Beta(1, 1)
+              def transition = () => {
+                val next = Universe.createNew
+                val f = Binomial(5,b)("binomial",next)
+                next
+              }
+              val algorithm = EMWithBP.online(transition, b)(initial)
+              algorithm.start()
+              algorithm.update(List(NamedEvidence("binomial", Observation(4))))
+              algorithm.update(List(NamedEvidence("binomial", Observation(3))))
+              val result = b.MAPValue
+              algorithm.kill
+              result should be(0.7 +- 0.01)
+          }
+
+        }
+
+      "used to estimate a Dirichlet parameter with two concentration parameters" should
+        {
+
+          "detect bias after a large enough number of trials" in
+            {
+              val initial = Universe.createNew
+              val d = Dirichlet(2, 2)
+              def transition = () => {
+                val next = Universe.createNew
+                val s = Select(d, true, false)("s",next)
+                next
+              }
+              val algorithm = EMWithBP.online(transition, d)(initial)
+              algorithm.start()
+              for (i <- 1 to 7) {
+                algorithm.update(List(NamedEvidence("s", Observation(true))))
+              }
+              for (i <- 1 to 3) {
+                algorithm.update(List(NamedEvidence("s", Observation(false))))
+              }
+              val result = d.MAPValue
+              algorithm.kill
+              result(0) should be(0.6666 +- 0.01)
+
+            }
+
+          "take the prior concentration parameters into account" in
+            {
+              val initial = Universe.createNew
+              val d = Dirichlet(3, 7)
+              def transition = () => {
+                val next = Universe.createNew
+                val s = Select(d, true, false)("s",next)
+                next
+              }
+              val algorithm = EMWithBP.online(transition, d)(initial)
+                         algorithm.start()
+              for (i <- 1 to 7) {
+                algorithm.update(List(NamedEvidence("s", Observation(true))))
+              }
+              for (i <- 1 to 3) {
+                algorithm.update(List(NamedEvidence("s", Observation(false))))
+              }
+              val result = d.MAPValue
+              algorithm.kill
+              result(0) should be(0.50 +- 0.01)
+
+            }
+
+        }        
+          
+          "used to estimate multiple parameters" should
+            {
+
+              "leave parameters having no observations unchanged" in
+                {
+                  val initial = Universe.createNew
+                  val d = Dirichlet(2.0, 4.0, 2.0)
+                  val b = Beta(2.0, 2.0)
+                  val outcomes = List(1, 2, 3)
+
+                  def transition = () => {
+                    val next = Universe.createNew
+                    val s = Select(d, outcomes:_*)("s",next)
+                    next
+                  }
+                  val algorithm = EMWithBP.online(transition,d,b)(initial)
+                             algorithm.start()
+                  for (i <- 1 to 4) {
+                    algorithm.update(List(NamedEvidence("s", Observation(1))))
+                  }
+
+                  for (i <- 1 to 2) {
+                    algorithm.update(List(NamedEvidence("s", Observation(2))))
+                  }
+
+                  for (i <- 1 to 4) {
+                    algorithm.update(List(NamedEvidence("s", Observation(3))))
+                  }
+
+               
+
+                  val result = d.MAPValue
+         
+                  result(0) should be(0.33 +- 0.01)
+                  result(1) should be(0.33 +- 0.01)
+                  result(2) should be(0.33 +- 0.01)
+
+                  val betaResult = b.MAPValue
+                  betaResult should be(0.5)
+                  algorithm.kill
+                }
+
+              "correctly estimate all parameters with observations" in
+                {
+                  val initial = Universe.createNew
+                  val d = Dirichlet(2.0, 3.0, 2.0)
+                  val b = Beta(3.0, 7.0)
+                  val outcomes = List(1, 2, 3)
+                  
+                  def transition = () => {
+                    val next = Universe.createNew
+                    val s = Select(d, outcomes:_*)("s",next)
+                    val f = Flip(b)("f",next)
+                    next
+                  }
+                  val algorithm = EMWithBP.online(transition,d,b)(initial)
+                             algorithm.start()
+                  for (i <- 1 to 4) {
+                    algorithm.update(List(NamedEvidence("s", Observation(1))))
+                  }
+
+                  for (i <- 1 to 2) {
+                    algorithm.update(List(NamedEvidence("s", Observation(2))))
+                  }
+
+                  for (i <- 1 to 4) {
+                    algorithm.update(List(NamedEvidence("s", Observation(3))))
+                  }
+                  
+                  for (i <- 1 to 7) {
+                    algorithm.update(List(NamedEvidence("f", Observation(true))))
+                  }
+                  
+                  for (i <- 1 to 3) {
+                    algorithm.update(List(NamedEvidence("f", Observation(false))))
+                  }
+
+                  val result = d.MAPValue
+
+                  result(0) should be(0.35 +- 0.01)
+                  result(1) should be(0.28 +- 0.01)
+                  result(2) should be(0.36 +- 0.01)
+
+                  val betaResult = b.MAPValue
+                  betaResult should be(0.5 +- 0.01)
+                 algorithm.kill
+                }
+            }
+    }
+  
   "Online expectation maximization" when
     {
 
diff --git a/Figaro/src/test/scala/com/cra/figaro/test/algorithm/sampling/AnnealingTest.scala b/Figaro/src/test/scala/com/cra/figaro/test/algorithm/sampling/AnnealingTest.scala
index 83691621..4a7cf5e3 100644
--- a/Figaro/src/test/scala/com/cra/figaro/test/algorithm/sampling/AnnealingTest.scala
+++ b/Figaro/src/test/scala/com/cra/figaro/test/algorithm/sampling/AnnealingTest.scala
@@ -22,6 +22,9 @@ import com.cra.figaro.language._
 import com.cra.figaro.library.atomic.continuous._
 import com.cra.figaro.library.compound.^^
 import com.cra.figaro.test._
+import scala.util.Random
+import com.cra.figaro.library.atomic.discrete.SwitchingFlip
+import com.cra.figaro.library.compound.If
 
 class AnnealingTest extends WordSpec with Matchers with PrivateMethodTester {
 
@@ -63,7 +66,7 @@ class AnnealingTest extends WordSpec with Matchers with PrivateMethodTester {
     "converge to the most likely value with an interval using anytime" in {
       Universe.createNew()
       val elems = buildUndirected(4, Flip(0.5), (b: (Boolean, Boolean)) => if (b._1 && b._2) 2.0 else 1.0)
-      val annealer = MetropolisHastingsAnnealer(ProposalScheme.default, Schedule.default(2.0),1,2)
+      val annealer = MetropolisHastingsAnnealer(ProposalScheme.default, Schedule.default(2.0), 1, 2)
       annealer.start()
       Thread.sleep(500)
       annealer.stop()
@@ -73,11 +76,11 @@ class AnnealingTest extends WordSpec with Matchers with PrivateMethodTester {
       }
       annealer.kill
     }
-    
+
     "converge to the most likely value with an interval using one-time" in {
       Universe.createNew()
       val elems = buildUndirected(4, Flip(0.5), (b: (Boolean, Boolean)) => if (b._1 && b._2) 2.0 else 1.0)
-      val annealer = MetropolisHastingsAnnealer(10000, ProposalScheme.default, Schedule.default(2.0),1,2)
+      val annealer = MetropolisHastingsAnnealer(10000, ProposalScheme.default, Schedule.default(2.0), 1, 2)
       annealer.start()
       annealer.stop()
       for { i <- 1 to 4 } {
@@ -86,40 +89,39 @@ class AnnealingTest extends WordSpec with Matchers with PrivateMethodTester {
       }
       annealer.kill
     }
-    
-    
+
     "produce higher temperatue with higher k" in {
       Universe.createNew()
       val elems1 = buildUndirected(4, Flip(0.5), (b: (Boolean, Boolean)) => if (b._1 && b._2) 2.0 else 1.0)
-      val annealer1 = MetropolisHastingsAnnealer(1000, ProposalScheme.default, Schedule.default(4.0),100)
+      val annealer1 = MetropolisHastingsAnnealer(1000, ProposalScheme.default, Schedule.default(4.0), 100)
       annealer1.start()
       val temp1 = annealer1.getTemperature
 
       val elems2 = buildUndirected(4, Flip(0.5), (b: (Boolean, Boolean)) => if (b._1 && b._2) 2.0 else 1.0)
-      val annealer2 = MetropolisHastingsAnnealer(1000, ProposalScheme.default, Schedule.default(2.0),100)
+      val annealer2 = MetropolisHastingsAnnealer(1000, ProposalScheme.default, Schedule.default(2.0), 100)
       annealer2.start()
       val temp2 = annealer2.getTemperature
 
       temp2 should be > temp1
-      
+
       annealer1.kill
       annealer2.kill
     }
-    
+
     "produce higher temperatue with higher k with burn-in" in {
       Universe.createNew()
       val elems1 = buildUndirected(4, Flip(0.5), (b: (Boolean, Boolean)) => if (b._1 && b._2) 2.0 else 1.0)
-      val annealer1 = MetropolisHastingsAnnealer(1000, ProposalScheme.default, Schedule.default(4.0),100)
+      val annealer1 = MetropolisHastingsAnnealer(1000, ProposalScheme.default, Schedule.default(4.0), 100)
       annealer1.start()
       val temp1 = annealer1.getTemperature
 
       val elems2 = buildUndirected(4, Flip(0.5), (b: (Boolean, Boolean)) => if (b._1 && b._2) 2.0 else 1.0)
-      val annealer2 = MetropolisHastingsAnnealer(1000, ProposalScheme.default, Schedule.default(2.0),100)
+      val annealer2 = MetropolisHastingsAnnealer(1000, ProposalScheme.default, Schedule.default(2.0), 100)
       annealer2.start()
       val temp2 = annealer2.getTemperature
 
       temp2 should be > temp1
-      
+
       annealer1.kill
       annealer2.kill
     }
@@ -142,11 +144,11 @@ class AnnealingTest extends WordSpec with Matchers with PrivateMethodTester {
       temp2 should be > temp1
       annealer1.kill
     }
-    
+
     "increase the temperature with more iterations, including with burn-in" in {
       Universe.createNew()
       val elems1 = buildUndirected(4, Flip(0.5), (b: (Boolean, Boolean)) => if (b._1 && b._2) 2.0 else 1.0)
-      val annealer1 = MetropolisHastingsAnnealer(ProposalScheme.default, Schedule.default(2.0),100)
+      val annealer1 = MetropolisHastingsAnnealer(ProposalScheme.default, Schedule.default(2.0), 100)
       annealer1.start()
       Thread.sleep(250)
       annealer1.stop
@@ -164,4 +166,39 @@ class AnnealingTest extends WordSpec with Matchers with PrivateMethodTester {
 
   }
 
+  "Running anytime annealing" should {
+
+    def buildComplicatedModel(): Universe = {
+      val u = Universe.createNew()
+      val flips = Seq.fill(1000)(Math.pow(0.01, Random.nextDouble)).map(SwitchingFlip(_))
+      val ifs = for ((f, i) <- flips.zipWithIndex) yield {
+        val fi = If(f, SwitchingFlip(Math.pow(0.01, Random.nextDouble)), Constant(false))(i.toString(), u)
+        if (Random.nextDouble > 0.80) fi.observe(Random.nextBoolean())
+      }
+      u
+    }
+
+    "give a most likely value at any time" in {
+      try {
+        val u = buildComplicatedModel()
+        val annealer = MetropolisHastingsAnnealer(ProposalScheme.default, Schedule.default(2.0))(u)
+        if (annealer.isActive) annealer.resume else annealer.start
+        Thread.sleep(100)
+        annealer.stop()
+        for (i <- 1 to 1000) {
+          if (annealer.isActive) annealer.resume else annealer.start
+          Thread.sleep(25)
+          annealer.stop()
+          for (i <- 0 until 1000) {
+            annealer.mostLikelyValue(u.getElementByReference[Boolean](i.toString()))
+          }
+        }
+        annealer.kill
+      } catch {
+        case knf: Exception => fail("running algorithm should not produce exceptions.")
+        case _ =>
+      }
+    }
+  }
+
 }
diff --git a/Figaro/src/test/scala/com/cra/figaro/test/algorithm/sampling/MHTest.scala b/Figaro/src/test/scala/com/cra/figaro/test/algorithm/sampling/MHTest.scala
index bab3bc92..28a6815d 100644
--- a/Figaro/src/test/scala/com/cra/figaro/test/algorithm/sampling/MHTest.scala
+++ b/Figaro/src/test/scala/com/cra/figaro/test/algorithm/sampling/MHTest.scala
@@ -32,7 +32,7 @@ import com.cra.figaro.test.tags.Performance
 
 class MHTest extends WordSpec with Matchers with PrivateMethodTester {
   "Anytime MetropolisHastings" should {
-        
+
     "for a constant produce the constant with probability 1" in {
       Universe.createNew()
       val c = Constant(1)
@@ -263,7 +263,7 @@ class MHTest extends WordSpec with Matchers with PrivateMethodTester {
       }
     }
 
-    "with a condition on the result of an If in which both consequences are random, correctly condition the test" taggedAs(NonDeterministic) in {
+    "with a condition on the result of an If in which both consequences are random, correctly condition the test" taggedAs (NonDeterministic) in {
       val numSamples = 100000
       val tolerance = 0.01
       Universe.createNew()
@@ -336,8 +336,8 @@ class MHTest extends WordSpec with Matchers with PrivateMethodTester {
   }
 
   "Testing a Metropolis-Hastings algorithm" should {
-    
-   "produce correct results with a typed scheme" in {
+
+    "produce correct results with a typed scheme" in {
       val numSamples = 200000
       val tolerance = 0.01
       Universe.createNew()
@@ -366,23 +366,22 @@ class MHTest extends WordSpec with Matchers with PrivateMethodTester {
         mh.kill()
       }
     }
-   
-   //Need  better use cases for these two schemes.
+
+    //Need  better use cases for these two schemes.
     "produce correct results with a switching scheme" in {
       Universe.createNew()
       val s1 = Select(0.4 -> 1, 0.6 -> 4)
       val s2 = Select(0.2 -> 2, 0.8 -> 3)
       s1.setConstraint((i: Int) => if (i == 1) 1; else 0.25)
       val p1 = Predicate(s1, (i: Int) => i == 4)
-      val scheme = SwitchScheme(() => (s1,s2),Some(FinalScheme(() => Universe.universe.randomStochasticElement())))
+      val scheme = SwitchScheme(() => (s1, s2), Some(FinalScheme(() => Universe.universe.randomStochasticElement())))
       val alg = MetropolisHastings(scheme)
       s1.value = 1
       val (_, scores, _) = alg.test(10000, List(p1), List())
       // probability of s1 changing to 4 is 0.5 (prob. s1 proposed) * 0.6 (prob 4 is chosen if s1 proposed) * 0.25
       scores(p1) should be(0.075 +- 0.01)
     }
-    
-    
+
     "produce values satisfying a predicate the correct amount of the time without conditions or constraints" in {
       Universe.createNew()
       val s1 = Select(0.4 -> 1, 0.6 -> 4)
@@ -421,7 +420,7 @@ class MHTest extends WordSpec with Matchers with PrivateMethodTester {
   }
 
   "Fixed bugs" should {
-    "propose temporary elements when non-temporary elements are non-stochastic" taggedAs(NonDeterministic) in {
+    "propose temporary elements when non-temporary elements are non-stochastic" taggedAs (NonDeterministic) in {
       Universe.createNew()
       val c1 = NonCachingChain(Constant(0), (i: Int) => Flip(0.7))
       val a1 = If(c1, Constant(1), Constant(0))
@@ -430,6 +429,47 @@ class MHTest extends WordSpec with Matchers with PrivateMethodTester {
       alg.start()
       alg.probability(a1, 1) should be(0.7 +- 0.01)
     }
+
+    "not record invalid states" taggedAs (NonDeterministic) in {
+      Universe.createNew()
+      val c1 = Flip(0.9)
+      val a1 = If(c1, com.cra.figaro.library.atomic.discrete.Uniform(1, 2), Select(0.999 -> 2, 0.001 -> 3))
+      a1.observe(3)
+      val alg = MetropolisHastings(10000, ProposalScheme.default, a1, c1)
+      alg.start()
+      alg.distribution(c1).toList.exists(s => s._2 == true) should be(false)
+    }
+    "correctly update elements in order" taggedAs (NonDeterministic) in {
+      // bug test for issue #453
+
+      Universe.createNew()
+      val Market = Select(0.5 -> 0, 0.3 -> 1, 0.2 -> 2)
+      val Survey = CPD(Constant(true), Market,
+        (false, 0) -> Constant(-1), (false, 1) -> Constant(-1), (false, 2) -> Constant(-1),
+        (true, 0) -> Select(0.6 -> 0, 0.3 -> 1, 0.1 -> 2),
+        (true, 1) -> Select(0.3 -> 0, 0.4 -> 1, 0.3 -> 2),
+        (true, 2) -> Select(0.1 -> 0, 0.4 -> 1, 0.5 -> 2))
+      val Found = Chain(Survey, (i: Int) => com.cra.figaro.library.atomic.discrete.Uniform(true, false))
+      def Value_fcn(f: Boolean, m: Int): Double = {
+        if (f) {
+          m match {
+            case 0 => -7.0
+            case 1 => 5.0
+            case 2 => 20.0
+          }
+        } else {
+          0.0
+        }
+      }
+      val Value = Apply(Found, Market, Value_fcn)
+      val Z = ^^(Found, Value)
+
+      val alg = MetropolisHastings(20000, ProposalScheme.default, Z)
+      alg.start()
+
+      val a = alg.distribution(Z).toList
+      alg.distribution(Z).toList.exists(s => s._2._1 == false && s._2._2 != 0.0) should be(false)
+    }
   }
 
   def newState(mh: MetropolisHastings): State = {
@@ -508,7 +548,7 @@ class MHTest extends WordSpec with Matchers with PrivateMethodTester {
       mh.start()
       Thread.sleep(1000)
       mh.stop()
-//      println(mh.getSampleCount.toString + " samples")
+      //      println(mh.getSampleCount.toString + " samples")
       mh.probability(target, predicate) should be(prob +- tolerance)
     } finally {
       mh.kill()
diff --git a/Figaro/src/test/scala/com/cra/figaro/test/example/MultiDecisionTest.scala b/Figaro/src/test/scala/com/cra/figaro/test/example/MultiDecisionTest.scala
index 72fe3e06..61953081 100644
--- a/Figaro/src/test/scala/com/cra/figaro/test/example/MultiDecisionTest.scala
+++ b/Figaro/src/test/scala/com/cra/figaro/test/example/MultiDecisionTest.scala
@@ -58,13 +58,13 @@ class MultiDecisionTest extends WordSpec with Matchers {
         override def oneTest = {
 
           val result = doTest((e1: List[Element[Double]], e2: List[Decision[_, _]]) =>
-            MultiDecisionMetropolisHastings(200000, propmaker, 10000, e1, e2: _*))
+            MultiDecisionMetropolisHastings(300000, propmaker, 20000, e1, e2: _*))
 
           update(result(0), NDTest.BOOLEAN, "MHMulti-DecisionFound(-1)", true, .90)
           update(result(1), NDTest.BOOLEAN, "MHMulti-DecisionFound(0)", false, .90)
           update(result(2), NDTest.BOOLEAN, "MHMulti-DecisionFound(1)", true, .90)
           update(result(3), NDTest.BOOLEAN, "MHMulti-DecisionFound(2)", true, .90)
-          update(result(4), NDTest.BOOLEAN, "MHMulti-DecisionTest(0)", true, .90)
+          update(result(4), NDTest.BOOLEAN, "MHMulti-DecisionTest(0)", true, .70)
         }
       }
       ndtest.run(10)
diff --git a/Figaro/src/test/scala/com/cra/figaro/test/example/MultiValuedTest.scala b/Figaro/src/test/scala/com/cra/figaro/test/example/MultiValuedTest.scala
index 34e3cf8b..144a7ece 100644
--- a/Figaro/src/test/scala/com/cra/figaro/test/example/MultiValuedTest.scala
+++ b/Figaro/src/test/scala/com/cra/figaro/test/example/MultiValuedTest.scala
@@ -1,13 +1,13 @@
 /*
- * MultiValuedTest.scala 
+ * MultiValuedTest.scala
  * Multi-valued reference uncertainty example tests.
- * 
+ *
  * Created By:      Avi Pfeffer (apfeffer@cra.com)
  * Creation Date:   Jan 1, 2009
- * 
+ *
  * Copyright 2013 Avrom J. Pfeffer and Charles River Analytics, Inc.
  * See http://www.cra.com or email figaro@cra.com for information.
- * 
+ *
  * See http://www.github.com/p2t2/figaro for a copy of the software license.
  */
 
diff --git a/Figaro/src/test/scala/com/cra/figaro/test/example/MutableMovieTest.scala b/Figaro/src/test/scala/com/cra/figaro/test/example/MutableMovieTest.scala
index f781b137..4d330d43 100644
--- a/Figaro/src/test/scala/com/cra/figaro/test/example/MutableMovieTest.scala
+++ b/Figaro/src/test/scala/com/cra/figaro/test/example/MutableMovieTest.scala
@@ -35,8 +35,7 @@ class MutableMovieTest extends WordSpec with Matchers {
     }
 
     "produce the correct probability under Metropolis-Hastings" taggedAs (Example, NonDeterministic) in {
-      test((e: Element[Boolean]) => MetropolisHastings(200000, chooseScheme, 10000, e))
-      println(com.cra.figaro.util.seed)
+      test((e: Element[Boolean]) => MetropolisHastings(300000, chooseScheme, 20000, e))      
     }
   }
 
diff --git a/Figaro/src/test/scala/com/cra/figaro/test/experimental/structured/FactorMakerTest.scala b/Figaro/src/test/scala/com/cra/figaro/test/experimental/structured/FactorMakerTest.scala
index 1f857968..9b96bcb2 100644
--- a/Figaro/src/test/scala/com/cra/figaro/test/experimental/structured/FactorMakerTest.scala
+++ b/Figaro/src/test/scala/com/cra/figaro/test/experimental/structured/FactorMakerTest.scala
@@ -3588,7 +3588,7 @@ class FactorMakerTest extends WordSpec with Matchers {
         factor1L.get(List(e2Index2, 0)) should equal (1.0)
         factor2L.size should equal (2)
         factor2L.get(List(e2IndexStar, 0)) should equal (0.3 +- .0001)
-        factor2L.get(List(e2Index2, 0)) should equal (0.3)
+        factor2L.get(List(e2Index2, 0)) should equal (0.3 +- .0001)
         val List(factor1U) = c11.constraintUpper
         val List(factor2U) = c12.constraintUpper
         factor1U.variables should equal (List(c2.variable, c11.variable))
@@ -3598,7 +3598,7 @@ class FactorMakerTest extends WordSpec with Matchers {
         factor1U.get(List(e2Index2, 0)) should equal (1.0)
         factor2U.size should equal (2)
         factor2U.get(List(e2IndexStar, 0)) should equal (1.0)
-        factor2U.get(List(e2Index2, 0)) should equal (0.3)
+        factor2U.get(List(e2Index2, 0)) should equal (0.3 +- .0001)
       }
     }
 
diff --git a/FigaroExamples/META-INF/MANIFEST.MF b/FigaroExamples/META-INF/MANIFEST.MF
index 9c4a0ab6..45a4a81e 100644
--- a/FigaroExamples/META-INF/MANIFEST.MF
+++ b/FigaroExamples/META-INF/MANIFEST.MF
@@ -2,8 +2,8 @@ Manifest-Version: 1.0
 Bundle-ManifestVersion: 2
 Bundle-Name: FigaroExamples
 Bundle-SymbolicName: com.cra.figaro.examples
-Bundle-Version: 3.1.0
-Require-Bundle: com.cra.figaro;bundle-version="3.1.0",
+Bundle-Version: 3.2.0
+Require-Bundle: com.cra.figaro;bundle-version="3.2.0",
  org.scala-lang.scala-library
 Bundle-Vendor: Charles River Analytics
 Bundle-RequiredExecutionEnvironment: JavaSE-1.6
diff --git a/README.md b/README.md
index fc177931..383cab37 100644
--- a/README.md
+++ b/README.md
@@ -4,4 +4,4 @@ Figaro is a probabilistic programming language that supports development of very
 
 Figaro makes it possible to express probabilistic models using the power of programming languages, giving the modeler the expressive tools to create a wide variety of models. Figaro comes with a number of built-in reasoning algorithms that can be applied automatically to new models. In addition, Figaro models are data structures in the Scala programming language, which is interoperable with Java, and can be constructed, manipulated, and used directly within any Scala or Java program.
 
-Figaro is free and is released under an [open-source license](https://github.com/p2t2/figaro/blob/master/LICENSE). The current, stable binary release of Figaro can be found [here](https://www.cra.com/figaro). For more information please see the [Figaro Guide](https://www.cra.com/pdf/Figaro-Quick-Start-Guide.pdf) and [Figaro Tutorial](https://www.cra.com/pdf/Figaro-tutorial.pdf). Documentation of the Figaro library interface can be found [here](https://www.cra.com/files/Figaro_Scaladoc/index.html).
+Figaro is free and is released under an [open-source license](https://github.com/p2t2/figaro/blob/master/LICENSE). The current, stable binary release of Figaro can be found [here](https://www.cra.com/figaro). For more information please see the [Figaro Release Notes](https://www.cra.com/pdf/Figaro-Release-Notes.pdf) and [Figaro Tutorial](https://www.cra.com/pdf/Figaro-tutorial.pdf). Documentation of the Figaro library interface can be found [here](https://www.cra.com/files/Figaro_Scaladoc/index.html).
diff --git a/doc/Figaro Release Notes.pdf b/doc/Figaro Release Notes.pdf
index 2704c0e3..d7394855 100644
Binary files a/doc/Figaro Release Notes.pdf and b/doc/Figaro Release Notes.pdf differ
diff --git a/project/Build.scala b/project/Build.scala
index 26d3a612..28640f36 100644
--- a/project/Build.scala
+++ b/project/Build.scala
@@ -21,24 +21,11 @@ import com.typesafe.sbteclipse.plugin.EclipsePlugin.EclipseKeys
 
 object FigaroBuild extends Build {
 
-  // Copy dependency JARs to /target/<scala-version>/lib
-  // Courtesy of
-  // http://stackoverflow.com/questions/7351280/collecting-dependencies-under-sbt-0-10-putting-all-dependency-jars-to-target-sc
-  lazy val copyDependencies = TaskKey[Unit]("copy-deps")
-
-  def copyDepTask = copyDependencies <<= (update, crossTarget, scalaVersion) map {
-    (updateReport, out, scalaVer) =>
-    updateReport.allFiles foreach { srcPath =>
-      val destPath = out / "lib" / srcPath.getName
-      IO.copyFile(srcPath, destPath, preserveLastModified=true)
-    }
-  }
-
   override val settings = super.settings ++ Seq(
     organization := "com.cra.figaro",
     description := "Figaro: a language for probablistic programming",
-    version := "3.1.0.0",
-    scalaVersion := "2.11.4",
+    version := "3.2.0.0",
+    scalaVersion := "2.11.6",
     crossPaths := true,
     publishMavenStyle := true,
     pomExtra :=
@@ -101,6 +88,8 @@ object FigaroBuild extends Build {
       "com.storm-enroute" %% "scalameter" % "0.6" % "provided",
       "org.scalatest" %% "scalatest" % "2.2.4" % "provided, test"
     ))
+    // Copy all managed dependencies to \lib_managed directory
+    .settings(retrieveManaged := true)
     // Enable forking
     .settings(fork := true)
     // Increase max memory for JVM for both testing and runtime
@@ -122,8 +111,6 @@ object FigaroBuild extends Build {
 	val cp = (fullClasspath in assembly).value
 	cp filter {_.data.getName == "arpack_combined_all-0.1-javadoc.jar"}
     })
-    // Copy dependency JARs
-    .settings(copyDepTask)
     // ScalaMeter settings
     .settings(testFrameworks += new TestFramework("org.scalameter.ScalaMeterFramework"))
     .settings(logBuffered := false)
@@ -133,8 +120,6 @@ object FigaroBuild extends Build {
   lazy val examples = Project("FigaroExamples", file("FigaroExamples"))
     .dependsOn(figaro)
     .settings(packageOptions := Seq(Package.JarManifest(examplesManifest)))
-    // Copy dependency JARs
-    .settings(copyDepTask)
     // SBTEclipse settings
     .settings(EclipseKeys.eclipseOutput := Some("target/scala-2.11/classes"))