Pull Request

.travis.yml

@@ -0,0 +1,11 @@
+language: python
+python:
+  - "2.7"
+  - "3.4"
+  - "3.5"
+  - "3.6"
+
+# command to install dependencies
+#install: "pip install -r requirements.txt"
+# command to run tests
+script: py.test

line.py

@@ -0,0 +1,153 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+from __future__ import print_function
+
+import emcee
+import corner
+import numpy as np
+import scipy.optimize as op
+import matplotlib.pyplot as pl
+from matplotlib.ticker import MaxNLocator
+
+# Reproducible results!
+np.random.seed(123)
+
+# Choose the "true" parameters.
+m_true = -0.9594
+b_true = 4.294
+f_true = 0.534
+
+# Generate some synthetic data from the model.
+
+y_list = []
+x_list = []
+yerr_list = []
+
+N = 50
+x = np.sort(10*np.random.rand(N))
+yerr = 0.1+0.5*np.random.rand(N)
+y = m_true*x+b_true
+y += np.abs(f_true*y) * np.random.randn(N)
+y += yerr * np.random.randn(N)
+
+# Plot the dataset and the true model.
+xl = np.array([0, 10])
+pl.errorbar(x, y, yerr=yerr, fmt=".k")
+pl.plot(xl, m_true*xl+b_true, "k", lw=3, alpha=0.6)
+pl.ylim(-9, 9)
+pl.xlabel("$x$")
+pl.ylabel("$y$")
+pl.tight_layout()
+pl.savefig("line-data.png")
+
+# Do the least-squares fit and compute the uncertainties.
+A = np.vstack((np.ones_like(x), x)).T
+C = np.diag(yerr * yerr)
+cov = np.linalg.inv(np.dot(A.T, np.linalg.solve(C, A)))
+b_ls, m_ls = np.dot(cov, np.dot(A.T, np.linalg.solve(C, y)))
+print("""Least-squares results:
+    m = {0} ± {1} (truth: {2})
+    b = {3} ± {4} (truth: {5})
+""".format(m_ls, np.sqrt(cov[1, 1]), m_true, b_ls, np.sqrt(cov[0, 0]), b_true))
+
+# Plot the least-squares result.
+print(m_ls*xl+b_ls,xl,"this is it")
+pl.plot(xl, m_ls*xl+b_ls, "--k")
+pl.savefig("line-least-squares.png")
+
+# Define the probability function as likelihood * prior.
+def lnprior(theta):
+    m, b, lnf = theta
+    if -5.0 < m < 0.5 and 0.0 < b < 10.0 and -10.0 < lnf < 1.0:
+        return 0.0
+    return -np.inf
+
+def lnlike(theta, x, y, yerr):
+    m, b, lnf = theta
+    model = m * x + b
+    inv_sigma2 = 1.0/(yerr**2 + model**2*np.exp(2*lnf))
+    return -0.5*(np.sum((y-model)**2*inv_sigma2 - np.log(inv_sigma2)))
+
+def lnprob(theta, x, y, yerr):
+    lp = lnprior(theta)
+    if not np.isfinite(lp):
+        return -np.inf
+        print(lp + loglikechain[0])
+    return lp + lnlike(theta, x, y, yerr)
+
+# Find the maximum likelihood value.
+chi2 = lambda *args: -2 * lnlike(*args)
+result = op.minimize(chi2, [m_true, b_true, np.log(f_true)], args=(x, y, yerr))
+m_ml, b_ml, lnf_ml = result["x"]
+print("""Maximum likelihood result:
+    m = {0} (truth: {1})
+    b = {2} (truth: {3})
+    f = {4} (truth: {5})
+""".format(m_ml, m_true, b_ml, b_true, np.exp(lnf_ml), f_true))
+
+# Plot the maximum likelihood result.
+pl.plot(xl, m_ml*xl+b_ml, "k", lw=2)
+pl.savefig("line-max-likelihood.png")
+
+# Set up the sampler.
+ndim, nwalkers = 3, 100
+pos = [result["x"] + 1e-4*np.random.randn(ndim) for i in range(nwalkers)]
+sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(x, y, yerr))
+
+# Clear and run the production chain.
+print("Running MCMC...")
+sampler.run_mcmc(pos, 500, rstate0=np.random.get_state())
+print("Done.")
+
+pl.clf()
+fig, axes = pl.subplots(3, 1, sharex=True, figsize=(8, 9))
+axes[0].plot(sampler.chain[:, :, 0].T, color="k", alpha=0.4)
+axes[0].yaxis.set_major_locator(MaxNLocator(5))
+axes[0].axhline(m_true, color="#888888", lw=2)
+axes[0].set_ylabel("$m$")
+
+axes[1].plot(sampler.chain[:, :, 1].T, color="k", alpha=0.4)
+axes[1].yaxis.set_major_locator(MaxNLocator(5))
+axes[1].axhline(b_true, color="#888888", lw=2)
+axes[1].set_ylabel("$b$")
+
+axes[2].plot(np.exp(sampler.chain[:, :, 2]).T, color="k", alpha=0.4)
+axes[2].yaxis.set_major_locator(MaxNLocator(5))
+axes[2].axhline(f_true, color="#888888", lw=2)
+axes[2].set_ylabel("$f$")
+axes[2].set_xlabel("step number")
+
+fig.tight_layout(h_pad=0.0)
+fig.savefig("line-time.png")
+
+# Make the triangle plot.
+burnin = 50
+samples = sampler.chain[:, burnin:, :].reshape((-1, ndim))
+
+fig = corner.corner(samples, labels=["$m$", "$b$", "$\ln\,f$"],
+                      truths=[m_true, b_true, np.log(f_true)])
+fig.savefig("line-triangle.png")
+
+# Plot some samples onto the data.
+pl.figure()
+for m, b, lnf in samples[np.random.randint(len(samples), size=100)]:
+    pl.plot(xl, m*xl+b, color="k", alpha=0.1)
+pl.plot(xl, m_true*xl+b_true, color="r", lw=2, alpha=0.8)
+pl.errorbar(x, y, yerr=yerr, fmt=".k")
+pl.ylim(-9, 9)
+pl.xlabel("$x$")
+pl.ylabel("$y$")
+pl.tight_layout()
+pl.savefig("line-mcmc.png")
+
+# Compute the quantiles.
+samples[:, 2] = np.exp(samples[:, 2])
+m_mcmc, b_mcmc, f_mcmc = map(lambda v: (v[1], v[2]-v[1], v[1]-v[0]),
+                             zip(*np.percentile(samples, [16, 50, 84],
+                                                axis=0)))
+print("""MCMC result:
+    m = {0[0]} +{0[1]} -{0[2]} (truth: {1})
+    b = {2[0]} +{2[1]} -{2[2]} (truth: {3})
+    f = {4[0]} +{4[1]} -{4[2]} (truth: {5})
+""".format(m_mcmc, m_true, b_mcmc, b_true, f_mcmc, f_true))

mcmc_test/Log_of_planck_function.py

@@ -0,0 +1,88 @@
+from __future__ import print_function
+
+import emcee
+import corner
+import numpy as np
+import scipy.optimize as op
+import matplotlib.pyplot as pl
+from matplotlib.ticker import MaxNLocator
+from scipy.optimize import curve_fit
+import PyAstronomy
+from astropy.io import ascii
+
+np.seterr(divide='ignore', invalid='ignore')
+
+h = 6.626070040*(10**(-34))  #J*s
+c = 299792458 #m/s
+k = 1.38064852*(10**(-23))  #m^2*kg*s^-2*K^-1
+
+a = 2*h*(c**2)
+b = (h*c)/k
+
+def my_own_Planck(x,T):
+    #x is the wavelength
+    #returns spectral radiance per sr
+    return((a/x**5)*(1/((np.exp(b/(T*x))-1))))
+
+h = 6.626e-34
+c = 3.0e+8
+k = 1.38e-23
+
+def planck(wav, T):
+    a = 2.0*h*c**2
+    b = h*c/(wav*k*T)
+    intensity = a/ ( (wav**5) * (np.exp(b) - 1.0) )
+    return intensity
+
+def model(microns,Teff,logfactor):
+  wavelength = microns*1.0e-6
+  flux=np.empty([len(wavelength)])
+  logflux=np.empty([len(wavelength)])
+  for i in range(len(wavelength)):
+    # Flux is the Planck function
+    flux[i] = ( (2.0*con.h*con.c**2)/(wavelength[i]**5) )/( np.exp( (con.h*con.c)/(con.k*Teff*wavelength[i]) ) - 1.0 )
+    # So logflux (which is what we want) is the log of this
+    logflux[i] = logfactor + np.log10(flux[i])
+  return logflux
+
+X = np.linspace(0,25,256,endpoint=True)
+
+results = ascii.read("results.dat")
+
+# Temperatures of other Brown Dwarfs used for comparison
+MCMC_Temp = results[0][0]
+T1, T1_name = 227, "WISE 1506+7027"
+T2, T2_name = 450, "WISE 0410+1502"
+T3, T3_name = 350, "WISE 1541−2250"
+
+# temperature and Brown Dwarf decriptions
+T1_description = T1_name + " (" + str(T1)+"K" + ")"
+T2_description = T1_name + " (" + str(T2)+"K" + ")"
+T3_description = T1_name + " (" + str(T3)+"K" + ")"
+MCMC_description = "WISE 0855-0714" + " (" + str(MCMC_Temp)+"K" + ")"
+
+# Ys used for plotting
+Y1 = (np.log10(planck(X*10**-6,T1)))
+Y2 = (np.log10(planck(X*10**-6,T2)))
+Y3 = (np.log10(planck(X*10**-6,T3)))
+Y4 = (np.log10(planck(X*10**-6,MCMC_Temp)))
+
+# Given data for WISE 0855-0714
+data = ascii.read("SED.dat", data_start=4)
+x = data[0][:]
+y = data[1][:]
+yerr = data[2][:]
+
+# Plotting the error bars for WISE 0855-0714
+pl.errorbar(x, y, yerr=yerr, fmt=".k")
+
+# Plotting the log of planck function for each brown dwarf
+plot1, = pl.plot(X,Y1, color="red", label=T1_description )
+plot2, = pl.plot(X,Y2, color="green", label=T2_description )
+plot3, = pl.plot(X,Y3, color="blue", label=T3_description )
+plot4, = pl.plot(X,Y4, color="black", label=MCMC_description )
+pl.legend([plot1,plot2,plot3,plot4],[T1_description, T2_description, T3_description, MCMC_description])
+pl.xlabel("Wavelength (um)")
+pl.ylabel("LOG10(Spectral Radiance)")
+
+pl.show()

mcmc_test/SED.dat

@@ -0,0 +1,11 @@
+Column #1 = Wavelength [microns]
+Column #2 = Log10( Flux [erg s^-1 cm^-2 Angstroms^-1] )
+Column #3 = Error on Column 2
+
+Col #1	Col #2	Col #3
+3.368	-18.435	0.087
+4.618	-17.042	0.087
+12.082	-15.728	0.087
+22.194	-16.307	0.087
+3.6	-18.063	0.043
+4.5	-17.173	0.043

mcmc_test/Spectral_Radiance.py

@@ -0,0 +1,71 @@
+from __future__ import print_function
+
+import emcee
+import corner
+import numpy as np
+import scipy.optimize as op
+import matplotlib.pyplot as pl
+from matplotlib.ticker import MaxNLocator
+from scipy.optimize import curve_fit
+import PyAstronomy
+import subprocess
+from astropy.io import ascii
+import os
+
+np.seterr(divide='ignore', invalid='ignore')
+np.seterr(over='ignore', invalid='ignore')
+
+results = ascii.read(os.path.dirname(os.path.realpath(__file__))+"\\results.dat")
+
+h = 6.626070040*(10**(-34))  #J*s
+c = 299792458 #m/s
+k_b = 1.38064852*(10**(-23))  #m^2*kg*s^-2*K^-1
+
+a = 2*h*(c**2)
+b = (h*c)/k_b
+
+def my_own_Planck(x,T):
+    #x is the wavelength
+    #returns spectral radiance per sr
+    return((a/x**5)*(1/((np.exp(b/(T*x))-1))))
+
+h = 6.626e-34
+c = 3.0e+8
+k = 1.38e-23
+
+def planck(wav, T):
+    a = 2.0*h*c**2
+    b = h*c/(wav*k*T)
+    intensity = a/ ( (wav**5) * (np.exp(b) - 1.0) )
+    return intensity
+
+X = np.linspace(0,50,1000,endpoint=True)
+
+# Temperatures of other Brown Dwarfs used for comparison
+MCMC_Temp = results[0][0]
+T1, T1_name = 227, "WISE 1506+7027"
+T2, T2_name = 450, "WISE 0410+1502"
+T3, T3_name = 350, "WISE 1541−2250"
+
+# temperature and Brown Dwarf decriptions
+T1_description = T1_name + " (" + str(T1)+"K" + ")"
+T2_description = T1_name + " (" + str(T2)+"K" + ")"
+T3_description = T1_name + " (" + str(T3)+"K" + ")"
+MCMC_description = "WISE 0855-0714" + " (" + str(MCMC_Temp)+"K" + ")"
+
+# Ys used for plotting
+Y1 = (planck(X*10**-6,T1))
+Y2 = (planck(X*10**-6,T2))
+Y3 = (planck(X*10**-6,T3))
+Y4 = (planck(X*10**-6,MCMC_Temp))
+
+# Plotting the planck function for each brown dwarf
+plot1, = pl.plot(X,Y1, color="red", label=T1_description )
+plot2, = pl.plot(X,Y2, color="green", label=T2_description )
+plot3, = pl.plot(X,Y3, color="blue", label=T3_description )
+plot4, = pl.plot(X,Y4, color="black", label=MCMC_description )
+pl.legend([plot1,plot2,plot3,plot4],[T1_description, T2_description, T3_description, MCMC_description])
+pl.xlabel("Wavelength (um)")
+pl.ylabel("Spectral Radiance")
+
+pl.show()