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>]

Research idea: NIR variation

I really want to know how near-infrared absorption lines correlate with emission lines; understanding this is essential for near infrared calibration.

The emission lines correlate reasonably well, but not perfectly.  There are models out there of how these should behave, but they seem to be proprietary and non-free so I have no interest in them - I won't pay before I know they work, and I can't test them.

But perhaps more careful measurements of the night sky lines could provide some information about the absorption, particularly in the Brackett-delta and Brackett-epsilon region we've been interested in lately.

I'm pretty sure the change in slope observed above is from observing at different airmasses (2.3 vs 1.05)

Research Idea: Stacking Finders

Idea: Stack all of the finders from spectroscopic observations.  Finder images tend to be on lower-quality CCDs with no filter, but they frequently produce very deep observations.  For example, the open K-band finder on TripleSpec (though it's technically not a CCD).

In order to stack them, you would need to mask out the bad pixels (already done) and compute astrometic solutions for the CCD.  Un-warping the images will take some work, but there should be plenty of information available from thousands of observations of different fields to make this computation nearly ideal.  Similarly, it should be possible to calibrate different pixels on the imager based on response to 2MASS standards.

Applications?  Very deep imaging of spectroscopic targets.  Short- and long-term variability (typical finder cadence is ~a few seconds).  Deep imaging around stars and galaxies of interest - probably far deeper than you could get with classical observing requests.

This project should be achievable by a motivated undergraduate, but I think the tools for astrometric solutions need to be in place first.  Astrometry.net is a great tool for this, but I think operates on spatial scales that are too large.  Once basic astrometric solutions are available (e.g., pointing center for the image), I think IRAF tools could be automated to compute the complete solution, which would then be applied to all images.  

Calibration might end up being the most challenging component, since there is variable atmospheric emission (absorption) that is not filtered by the finder.  Depending on the application, though, large calibration errors may be acceptable.  i.e., for deep nebular observations, morphology will be more important than absolute brightness, since the line responsible for the brightness cannot be directly determined.  Whereas, for variability, calibration is important, but it can be computed directly from other stars in the field.