@@ -538,12 +538,24 @@ template <class T, class MatA, class MatB>

 OptionalFloat GibbsSampler<T, MatA, MatB>::gibbsMass(AlphaParameters alpha,

 GapsRng *rng)

 {        

-    alpha *= mAnnealingTemp;

+    alpha.s *= mAnnealingTemp;

+    alpha.su *= mAnnealingTemp;

+

     if (alpha.s > gaps::epsilon)

     {

         float mean = (alpha.su - mLambda) / alpha.s;

         float sd = 1.f / std::sqrt(alpha.s);

-        return rng->truncNormal(gaps::epsilon, mMaxGibbsMass / mLambda, mean, sd);

+        float pLower = gaps::p_norm(0.f, mean, sd);

+

+        if (pLower < 1.f)

+        {

+            float m = rng->inverseNormSample(pLower, 1.f, mean, sd);

+            float gMass = gaps::min(m, mMaxGibbsMass / mLambda);

+            if (gMass >= gaps::epsilon)

+            {

+                return OptionalFloat(gMass);

+            }

+        }

     }

     return OptionalFloat();

 }

@@ -552,12 +564,22 @@ template <class T, class MatA, class MatB>

 OptionalFloat GibbsSampler<T, MatA, MatB>::gibbsMass(AlphaParameters alpha,

 float m1, float m2, GapsRng *rng)

 {

-    alpha *= mAnnealingTemp;

+    alpha.s *= mAnnealingTemp;

+    alpha.su *= mAnnealingTemp;

+

     if (alpha.s > gaps::epsilon)

     {

         float mean = alpha.su / alpha.s; // lambda cancels out

         float sd = 1.f / std::sqrt(alpha.s);

-        return rng->truncNormal(-m1, m2, mean, sd);

+        float pLower = gaps::p_norm(-m1, mean, sd);

+        float pUpper = gaps::p_norm(m2, mean, sd);

+

+        if (!(pLower >  0.95f || pUpper < 0.05f))

+        {

+            float delta = rng->inverseNormSample(pLower, pUpper, mean, sd);

+            float gMass = gaps::min(gaps::max(-m1, delta), m2); // conserve mass

+            return OptionalFloat(gMass);

+        }

     }

     return OptionalFloat();

 }

@@ -594,4 +616,4 @@ Archive& operator>>(Archive &ar, GibbsSampler<T, MatA, MatB> &s)

     return ar;

 }

-#endif // __COGAPS_GIBBS_SAMPLER_H__

+#endif // __COGAPS_GIBBS_SAMPLER_H__

\ No newline at end of file

@@ -227,23 +227,14 @@ float GapsRng::exponential(float lambda)

     return -1.f * std::log(uniform()) / lambda;

 }

-OptionalFloat GapsRng::truncNormal(float a, float b, float mean, float sd)

+float GapsRng::inverseNormSample(float a, float b, float mean, float sd)

 {

-    float pLower = gaps::p_norm(a, mean, sd);

-    float pUpper = gaps::p_norm(b, mean, sd);

-

-    if (!(pLower >  0.95f || pUpper < 0.05f))

+    float u = uniform(a, b);

+    while (u == 0.f || u == 1.f)

     {

-        float u = uniform(pLower, pUpper);

-        while (u == 0.f || u == 1.f)

-        {

-            u = uniform(pLower, pUpper);

-        }

-        float ret = gaps::q_norm(u, mean, sd);

-        ret = gaps::min(gaps::max(a, ret), b); // need this for truncation error

-        return OptionalFloat(ret);

+        u = uniform(a, b);

     }

-    OptionalFloat();

+    return gaps::q_norm(u, mean, sd);

 }

 float GapsRng::truncGammaUpper(float b, float shape, float scale)

@@ -319,4 +310,4 @@ float gaps::d_norm_fast(float d, float mean, float sd)

 {

     return std::exp((d - mean) * (d - mean) / (-2.f * sd * sd))

         / std::sqrt(2.f * gaps::pi * sd * sd);

-}

+}

\ No newline at end of file

@@ -72,7 +72,7 @@ public:

     int poisson(double lambda);

     float exponential(float lambda);

-    OptionalFloat truncNormal(float a, float b, float mean, float sd);

+    float inverseNormSample(float a, float b, float mean, float sd);

     float truncGammaUpper(float b, float shape, float scale);

     static void setSeed(uint64_t seed);

@@ -95,4 +95,4 @@ private:

     friend Archive& operator>>(Archive &ar, GapsRng &gen);

 };

-#endif // __COGAPS_RANDOM_H__

+#endif // __COGAPS_RANDOM_H__

\ No newline at end of file