To further elaborate my comment above as a potential answer to this question:
You are right in the sense that the NDF gives the statistical probability that the microsurface normal M exactly aligns with N (where exactly means: with respect to an infinitesimal solid angle around N and M).
As such, the NDF is a distribution function (more broadly a "generalized function") with unit one over steradians. That means, you get a unitless fraction value when you multiply the function's value with a the solid angle.
And this is where integration comes into play.
Because, what the NDF is ultimatly being used for is in the rendering equation to determine the fraction of light being reflected by the microfacets within the differential solid angle across M, where M is the half-vector between two other vectors V_o and V_i and where one of those two vectors is defined as the "outgoing direction" to a receiving differential surface and the other as the "incoming light direction".
We will use the w_o = outgoing and w_i = incoming direction, but due to reciprocity it really does not matter.
When we solve the rendering equation, we have to integrate over the whole hemisphere of the surface point we want to shade (i.e. whose radiance we want to compute along the outgoing direction w_o) to account for all light which will potentially be reflected off this surface point along the outgoing direction w_o. To do this, we have to know this point's "irradiance". And in order to get this, we evaluate the BRDF for all possible incoming light directions w_i and sum and average their contributions.
If we actually did sample all infinitely many w_i directions for a given w_o direction and asked the NDF for each such direction and averaged the NDF's cosine-weighted results, the result would always be exactly 1.0, because if we only consider the NDF as a potential factor in the BRDF, all outgoing light along w_o would have to come from somewhere when reflected off the surface.
And because of the Helmholtz reciprocity, we can rephrase this the other way around by saying: When fixing the w_i direction (so only considering a single direction light comes in from) and integrate over all possible outgoing directions w_o and again evaluated and averaged the cosine-weighted NDF, the result will yet again be exactly 1.0.
Because, when we again only consider the NDF in the BRDF, the incoming light will always be reflected to somewhere and not get lost.
So, after having writte all this, what does a single sample of the NDF function now actually mean? In short: It means you do not integrate over all differential solid angles and average the results, but only take a single sample for a single direction and you declare that as the result of your integral (which you only computed with a single sample). So you basically computed the Riemann Sum with a single sample.
Or in other words: Simply sampling the NDF at a single location carries no plausible physical meaning/result.
I hope that clears up some unknowns.
This was actually my first serious StackOverflow answer and I know I definitely can improve on writing those in the future!
Thanks! :)
