diff --git a/src/pygambit/nash.pxi b/src/pygambit/nash.pxi
index 1a681f8ffc..2e202146df 100644
--- a/src/pygambit/nash.pxi
+++ b/src/pygambit/nash.pxi
@@ -791,9 +791,7 @@ class LogitQREMixedBehaviorProfile:
     @property
     def game(self) -> Game:
         """The game on which this mixed strategy profile is defined."""
-        g = Game()
-        g.game = deref(self.thisptr).GetGame()
-        return g
+        return Game.wrap(deref(self.thisptr).GetGame())
 
     @property
     def lam(self) -> double:
@@ -808,11 +806,9 @@ class LogitQREMixedBehaviorProfile:
     @property
     def profile(self) -> MixedBehaviorProfileDouble:
         """The mixed strategy profile."""
-        profile = MixedBehaviorProfileDouble()
-        profile.profile = (
+        return MixedBehaviorProfileDouble.wrap(
             make_shared[c_MixedBehaviorProfile[double]](deref(self.thisptr).GetProfile())
         )
-        return profile
 
 
 def _logit_behavior_estimate(profile: MixedBehaviorProfileDouble,
diff --git a/src/pygambit/qre.py b/src/pygambit/qre.py
index b9377e5f14..0e2dd127c0 100644
--- a/src/pygambit/qre.py
+++ b/src/pygambit/qre.py
@@ -251,7 +251,9 @@ def _estimate_behavior_empirical(
     for player in data.game.players:
         for infoset in player.infosets:
             node = next(iter(infoset.members))
-            profile[node] = [math.exp(next(log_probs)) for _ in infoset.actions]
+            profile[node] = {
+                a.label: math.exp(next(log_probs)) for a in infoset.actions
+            }
     return LogitQREMixedBehaviorFitResult(
         data, "empirical", res.x[0], profile, -res.fun
     )
diff --git a/src/solvers/logit/efglogit.cc b/src/solvers/logit/efglogit.cc
index 1ee0b42d1e..ba9e1b5f1c 100644
--- a/src/solvers/logit/efglogit.cc
+++ b/src/solvers/logit/efglogit.cc
@@ -145,7 +145,7 @@ class EquationSystem {
   };
 
   Array<std::shared_ptr<Equation>> m_equations;
-  const Game &m_game;
+  Game m_game;
 };
 
 EquationSystem::EquationSystem(const Game &p_game) : m_game(p_game)
diff --git a/src/solvers/logit/nfglogit.cc b/src/solvers/logit/nfglogit.cc
index dc57aa3026..d838a2b449 100644
--- a/src/solvers/logit/nfglogit.cc
+++ b/src/solvers/logit/nfglogit.cc
@@ -220,7 +220,7 @@ class EquationSystem {
 
 private:
   std::vector<std::shared_ptr<Equation>> m_equations;
-  const Game &m_game;
+  Game m_game;
   mutable MixedStrategyProfile<double> m_profile, m_logProfile;
   mutable Vector<double> m_strategyValues;
   mutable Matrix<double> m_strategyDerivs;
diff --git a/tests/test_qre.py b/tests/test_qre.py
index 46d2e66353..419ada3a2a 100644
--- a/tests/test_qre.py
+++ b/tests/test_qre.py
@@ -2,11 +2,23 @@
 
 import pygambit as gbt
 
+from . import games
+
 
 def _asymmetric_2x2() -> gbt.Game:
     return gbt.Game.from_arrays([[1, 2], [3, 4]], [[4, 3], [2, 1]])
 
 
+def _asymmetric_poker_behavior_data() -> gbt.MixedBehaviorProfile:
+    game = games.create_stripped_down_poker_efg()
+    data = game.mixed_behavior_profile(rational=False)
+    for player in game.players:
+        for infoset in player.infosets:
+            node = next(iter(infoset.members))
+            data[node] = {a.label: float(i + 2) for i, a in enumerate(infoset.actions)}
+    return data
+
+
 @pytest.mark.parametrize("use_empirical", [False, True])
 def test_logit_estimate_strategy_accepts_rational_profile(use_empirical: bool):
     """logit_estimate() accepts a MixedStrategyProfileRational, converting it to
@@ -38,3 +50,31 @@ def test_logit_estimate_strategy_rational_and_float_data_agree():
                 rational_result.profile[player.label][strategy.label]
                 == pytest.approx(float_result.profile[player.label][strategy.label])
             )
+
+
+@pytest.mark.parametrize("use_empirical", [False, True])
+@pytest.mark.parametrize("local_max", [False, True])
+def test_logit_estimate_behavior_completes(use_empirical: bool, local_max: bool):
+    """logit_estimate() on a MixedBehaviorProfile returns a well-formed result instead of
+    crashing or raising.
+
+    Regression coverage for three bugs found while implementing rational-profile support
+    for `liap_solve`/`logit_estimate` (gambitproject/gambit#721, #458):
+    - `_estimate_behavior_fixedpoint()` segfaulted, because `EquationSystem` (in both
+      efglogit.cc and nfglogit.cc) stored its `Game` by dangling reference rather than
+      by value; `LogitBehaviorEstimate()` bound that reference straight to a temporary.
+    - `LogitQREMixedBehaviorProfile.profile`/`.game` (nash.pxi) constructed their return
+      value by calling the blocked `MixedBehaviorProfileDouble()`/`Game()` constructors
+      directly instead of via their `.wrap()` factories, raising `ValueError`.
+    - `_estimate_behavior_empirical()` assigned a `list` into a `MixedBehaviorProfile`,
+      which requires a `Mapping` from action label to weight, raising `TypeError`.
+    """
+    data = _asymmetric_poker_behavior_data()
+    result = gbt.qre.logit_estimate(data, use_empirical=use_empirical, local_max=local_max)
+    assert isinstance(result.profile, gbt.MixedBehaviorProfileDouble)
+    for player in data.game.players:
+        for infoset in player.infosets:
+            node = next(iter(infoset.members))
+            probs = dict(result.profile[node])
+            assert probs.keys() == {a.label for a in infoset.actions}
+            assert sum(probs.values()) == pytest.approx(1.0)