-->

Wednesday, September 19, 2012

Lognormal Density / IPython Notebook

An experiment with log normal density distributions (simple, pretty obvious) and a test of including ipython notebooks in blog posts. To do this, I created a gist (https://gist.github.com/3750781) and used the ipython notebook viewer (http://nbviewer.ipython.org/3750781/) to render the notebook in HTML. It worked pretty well; all I had to do was c&p the source and then change 2 lines of the style (background color -> black, foreground color -> white for the markup... maybe I should do it for the rest too).

Comparison of the derived power spectra and density PDFs for lognormal power distributions. Obviously (if you think about it), the power spectra are always flat for randomly distributed signal. The PDFs change significantly.

In [29]:
img = np.zeros([128,128])
norm = np.random.randn(128,128)
import agpy
lognorm = 10**(norm)*0.1 # peak amplitude ~ 0.1 Jy
lognorm_noisy = lognorm + np.random.randn(128,128)*0.03 # add noise at the 0.03 Jy level (like Bolocam)
clf();
subplot(121); imshow(log10(lognorm)); colorbar()
subplot(122); imshow(log10(lognorm_noisy)); colorbar();
In [30]:
rr,zz = agpy.psds.power_spectrum(lognorm)
rrn,zzn = agpy.psds.power_spectrum(lognorm_noisy)
In [31]:
clf()
subplot(121); loglog(rr,zz,label="Lognormal PSD"); loglog(rrn,zzn,label="Noisy Lognormal PSD")
subplot(122); hist(lognorm.flat,bins=logspace(-3,3)); gca().set_xscale('log')
hist(lognorm_noisy.flat,bins=logspace(-3,3),color='r',alpha=0.5,histtype='stepfilled');
figure();
semilogy(rr,zz,label="Lognormal PSD"); semilogy(rrn,zzn,label="Noisy Lognormal PSD")
Out [31]:
[<matplotlib.lines.Line2D at 0x111bd8390>]

So apparently there's no difference in the power spectrum and only a subtle difference in the PSD. I'll do the same thing with a lower peak, sub-noise

In [33]:
lognorm = 10**(norm)*0.01 # peak amplitude ~ 0.01 Jy
lognorm_noisy = lognorm + np.random.randn(128,128)*0.03 # add noise at the 0.03 Jy level (like Bolocam)
figure();
subplot(121); imshow(log10(lognorm)); colorbar()
subplot(122); imshow(log10(lognorm_noisy)); colorbar();
rr,zz = agpy.psds.power_spectrum(lognorm)
rrn,zzn = agpy.psds.power_spectrum(lognorm_noisy)
figure()
subplot(121); loglog(rr,zz,label="Lognormal PSD"); loglog(rrn,zzn,label="Noisy Lognormal PSD")
subplot(122); hist(lognorm.flat,bins=logspace(-3,3)); gca().set_xscale('log')
hist(lognorm_noisy.flat,bins=logspace(-3,3),color='r',alpha=0.5,histtype='stepfilled');
figure();
semilogy(rr,zz,label="Lognormal PSD"); semilogy(rrn,zzn,label="Noisy Lognormal PSD")
figure()
plot(rr,zz-zzn)
Out [33]:
[<matplotlib.lines.Line2D at 0x115609810>]

Any power spectrum differences are exceedingly subtle, at the $10^{-2}$ level

In [34]:
lognorm = 10**(norm)*0.001 # peak amplitude ~ 0.001 Jy
lognorm_noisy = lognorm + np.random.randn(128,128)*0.03 # add noise at the 0.03 Jy level (like Bolocam)
figure();
subplot(121); imshow(log10(lognorm)); colorbar()
subplot(122); imshow(log10(lognorm_noisy)); colorbar();
rr,zz = agpy.psds.power_spectrum(lognorm)
rrn,zzn = agpy.psds.power_spectrum(lognorm_noisy)
figure()
subplot(121); loglog(rr,zz,label="Lognormal PSD"); loglog(rrn,zzn,label="Noisy Lognormal PSD")
subplot(122); hist(lognorm.flat,bins=logspace(-3,3)); gca().set_xscale('log')
hist(lognorm_noisy.flat,bins=logspace(-3,3),color='r',alpha=0.5,histtype='stepfilled');
figure();
semilogy(rr,zz,label="Lognormal PSD"); semilogy(rrn,zzn,label="Noisy Lognormal PSD")
figure()
plot(rr,zz-zzn)
Out [34]:
[<matplotlib.lines.Line2D at 0x116ba1f10>]
In [35]:
lognorm = 10**(norm)*1.0 # peak amplitude ~ 1.0 Jy
lognorm_noisy = lognorm + np.random.randn(128,128)*0.03 # add noise at the 0.03 Jy level (like Bolocam)
figure();
subplot(121); imshow(log10(lognorm)); colorbar()
subplot(122); imshow(log10(lognorm_noisy)); colorbar();
rr,zz = agpy.psds.power_spectrum(lognorm)
rrn,zzn = agpy.psds.power_spectrum(lognorm_noisy)
figure()
subplot(121); loglog(rr,zz,label="Lognormal PSD"); loglog(rrn,zzn,label="Noisy Lognormal PSD")
subplot(122); hist(lognorm.flat,bins=logspace(-3,3)); gca().set_xscale('log')
hist(lognorm_noisy.flat,bins=logspace(-3,3),color='r',alpha=0.5,histtype='stepfilled');
figure();
semilogy(rr,zz,label="Lognormal PSD"); semilogy(rrn,zzn,label="Noisy Lognormal PSD")
figure()
plot(rr,zz-zzn)
Out [35]:
[<matplotlib.lines.Line2D at 0x118083f50>]

No comments:

Post a Comment