First, welcome to the comunity.
This is a tricky matter, as we don’t have access to the output vector… Well the output vector is easy to find in a glossy with no roughness, which is Out = 2*N*(N.In) - In (, or in 2.82+ series, with Out = -Reflect(I, N))… but with roughness, you don’t have N anymore; It is computed for each sample by microfaceting the original N.
To escape this dilema, you cannot use the builtin roughness; you need to create your own microfaceting, in order to have the real N (the half-vector to be more precise) for each sample. This however is not optimized for a clean rendering, and Cycles might miss a light if it’s too small, creating a huge amount of noise.
I use this method (explained here) to make diffraction gratings, thought it can be tweaked to produce other structural colors to some extend, but it requires lots of samples to get it look good.