Simply observe the noise for a longish while and estimate a PSD of it – for example, simply by doing an FFT and observing the magnitude of that, and calculating the mean square error to the theoretical (triangular) PSD of pink noise.
That can be easily implemented only using Python/numpy (fft, abs, mean are all implemented in numpy).
Another, pretty sensible, method would be designing a filter bank that exploits the fact that octaves always contain the same amount of power in pink noise (which is why you'd use it in e.g. audio room sounding).
scipy.signal
contains ample FIR filter design methods suitable for designing filters.
The band from 10-100 Hz should contain the same power (i.e. average sample square) as 100-1000 Hz, and 1000-10000 Hz, for example. Make sure scale the transition widths along with the passband widths. You don't have to go factors of 10, you can e.g. use 10-20 Hz, 20-40 Hz, 40-80 Hz … and verify that they all these contain the same power.