Mimicking Blender's "roughness" in my path tracer
For explicit light sampling: yep, you just evaluate the BRDF for that incoming direction and the output direction back toward the camera. There's probability involved in the case of an area light: you have to randomly choose a point on the light source, and include a factor to convert from probability density over light source area to a density over solid angle at the receiver point. It looks like you've got that already (that's the $\cos \theta' / (p(\mathbf{x}')\|\mathbf{x}-\mathbf{x}'\|^2)$ in the formula you mentioned).

For sampling the reflection ray, you could just pick a random ray from the hemisphere and weight it by the BRDF, the same as you'd do for diffuse. But that won't be very efficient; it'll take a huge number of samples to converge to a low-noise result.

What you really want to do is importance-sample the BRDF. This means choosing a random ray not uniformly over the hemisphere, but using a probability distribution that matches the shape of the BRDF (as much as possible). Then you'll automatically get rays that are clustered around the sharp reflection peak in the BRDF when the roughness is low, and spread out more over the hemisphere when the roughness is high.

So how do you do it? The paper Microfacet Models for Refraction through Rough Surfaces by Walter et al (which popularized the GGX distribution in graphics) handily gives importance-sampling formulas for several common BRDFs in section 5.2. Their approach is to importance-sample the normal distribution of the BRDF, then leave the Fresnel and geometry factors to be factored into the path weight.

For example, their GGX importance-sampling formula is (equations 35–36):
$$\begin{aligned}\theta_m &amp;= \arctan \left( \frac{\alpha_g\sqrt{\xi_1}}{\sqrt{1-\xi_1}} \right) \\ \phi_m &amp;= 2\pi\xi_2\end{aligned}$$
where $\alpha_g$ is the GGX roughness parameter, and $\xi_1, \xi_2$ are a pair of uniform random numbers in [0, 1]. The outputs $\theta_m, \phi_m$ are the polar coordinates (relative to the surface normal) of the sampled microfacet normal.

So to use this, you'd choose random inputs $\xi_1, \xi_2$, then evaluate the above equations, and use the resulting normal to compute the reflection ray, as well as the Fresnel and geometry factors from the BRDF.