--- a/src/modules/glm/samplers/HolmesHeld.h
+++ b/src/modules/glm/samplers/HolmesHeld.h
@@ -6,16 +6,68 @@
 namespace glm {
 
     /**
-     * Conjugate sampler for normal linear models.
+     * @short Holmes Held sampler for binary GLMs
+     *
+     * Sampler for probit and logistic regression models with binary
+     * outcome data, based on Holmes C and Held L (2006).  Bayesian
+     * Auxiliary Variables Models for Binary and Multinomial
+     * Regression, Bayesian Analysis, 1:148-168.
+     *
+     * In the Holmes-Held sampling method auxiliary normal variables
+     * (z[i]) are sampled from their joint distribution after
+     * marginalizing over the regression parameters of the GLM.  In
+     * probit models, this provides improved mixing over the
+     * Albert-Chib algorithm (See class AlbertChib). In logistic
+     * regression models, z[i] has a logistic distribution, which is
+     * represented by a scale mixture of normals.  The scale parameter
+     * (tau[i]) is sampled by the update method.  
+     *
+     * For logistic regression the improvement over the standard
+     * Albert-Chib algorithm is not so clear cut.
      */
     class HolmesHeld : public BinaryGLM {
-	bool _aux_init;
+	bool _aux_init; //Do we need to initialize auxiliary variables?
     public:
+	/**
+	 * Constructor.
+	 * 
+	 * @see GLMMethod#GLMMethod
+	 */
 	HolmesHeld(GraphView const *view, 
 		   std::vector<GraphView const *> const &sub_views,
 		   unsigned int chain);
+	/**
+	 * Updates the auxiliary variables (z[]) provided by the
+	 * parent BinaryGLM class. This function is called by
+	 * GLMMethod#updateLM.
+	 *
+	 * In the parent GLMMethod class, the posterior precision of
+	 * the regression parameters is represented by the matrix "A"
+	 * and the posterior mean "mu" solves A %*% mu = b.
+	 *
+	 * In this call, "N" holds the factorization 
+	 * P %*% A %*% t(P) = L %*% D %*% t(L)
+         * where P is a permutation matrix chosen to make
+	 * the factorization more efficient, and "w" solves
+	 * the equation L %*% w = P %*% b,
+	 *
+	 * These values are then used to calculate the marginal mean
+	 * and variance of z[], which forms a multivariate truncated
+	 * normal distribution. Each element of z[] is updated in turn
+	 * by Gibbs sampling.
+	 */
 	void updateAuxiliary(cholmod_dense *b, cholmod_factor *N, RNG *rng);
+	/**
+	 * Returns "Holmes-Held"
+	 */
 	std::string name() const;
+	/**
+	 * The update takes place in two steps. Firstly, for logistic
+	 * regression, the auxiliary variables tau[i] (representing
+	 * the precision of the latent variable z[i]) are updated
+	 * using the current values of z[], then GLMMethod#updateLM is
+	 * called. For probit regression, tau[i] is fixed at 1.
+	 */
 	void update(RNG *rng);
     };