<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Regression
This section follows the section 7.10 of Modern Statistical Methods in Astronomy by Feigelson and Babu
ordinary least squares
The statsmodels package provides an interface to fitting that is similar to what's available in R.
Step1:  The file SDSS_QSO.dat is from Feigelson's Astrostatistics school, September 2014 at Caltech
Step2:  Weighted least-squares
Step3:  Exercise
This example is drawn from Figure 8.2 of the AstroML book.
Generate data for 100 random supernovae.
Step4:  Put the data into a Pandas dataframe
Step5:  Fit a straight line
Step6:  Overplot the straight line fit
Step7:  In the cells below, fit 2nd and 4th degree polynomials and overplot on the data
<ASSISTANT_TASK:>
Python Code:
import statsmodels
statsmodels.__version__
import pandas as pd
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.sandbox.regression.predstd import wls_prediction_std
import matplotlib.pyplot as plt
Explanation: Regression
This section follows the section 7.10 of Modern Statistical Methods in Astronomy by Feigelson and Babu
ordinary least squares
The statsmodels package provides an interface to fitting that is similar to what's available in R.
End of explanation
df = pd.read_csv('data/SDSS_QSO.dat', delim_whitespace=True)
df.head()
df.describe()
# Keep only rows with positive magnitudes, and 18 < i < 22
df = df[(df.u_mag > 10) & (df.i_mag > 18) & (df.i_mag < 22) ]
# Set a limit on u magnitude errors
df.loc[df.sig_u_mag<0.02,'sig_u_mag'] = 0.02
df.sig_u_mag.min()
mod_ols = smf.ols('u_mag ~ i_mag', data=df)
res_ols = mod_ols.fit()
print(res_ols.summary())
print('Parameters: ', res_ols.params)
print('R2: ', res_ols.rsquared)
print('AIC', res_ols.aic)
Explanation: The file SDSS_QSO.dat is from Feigelson's Astrostatistics school, September 2014 at Caltech
End of explanation
mod_wls = smf.wls('u_mag ~ i_mag', data=df, weights=df.sig_u_mag**-2)
res_wls = mod_wls.fit()
print(res_wls.summary())
i_range = np.arange(18,22,0.2)
p_ols = res_ols.predict({"i_mag": i_range})
p_ols
%matplotlib inline
prstd, iv_l, iv_u = wls_prediction_std(res_wls,sm.add_constant(i_range), weights=1/(0.02*0.02))
fig, ax = plt.subplots(figsize=(8,6))
ax.scatter(df.i_mag, df.u_mag, marker='.', c='k', s=0.5, label="data", alpha=0.5)
ax.plot(i_range, res_wls.predict({"i_mag":i_range}), 'b-', label="Weighted LS")
ax.plot(i_range, res_ols.predict({"i_mag":i_range}), 'g--.', label="Ordinary LS")
ax.plot(i_range, iv_u, 'r--', label='95% confidence')
ax.plot(i_range, iv_l, 'r--')
ax.legend(loc='best')
plt.xlabel('SDSS i (mag)')
plt.ylabel('SDSS u (mag)')
plt.xlim(17.9, 21.6)
Explanation: Weighted least-squares
End of explanation
from astroML.cosmology import Cosmology
from astroML.datasets import generate_mu_z
import numpy as np
#------------------------------------------------------------
# Generate data
z_sample, mu_sample, dmu = generate_mu_z(100, random_state=0)
cosmo = Cosmology()
z = np.linspace(0.01, 2, 1000)
mu_true = np.asarray(map(cosmo.mu, z))
Explanation: Exercise
This example is drawn from Figure 8.2 of the AstroML book.
Generate data for 100 random supernovae.
End of explanation
sndata = pd.DataFrame({"z": z_sample, "mu": mu_sample, "mu_err": dmu})
fig = plt.figure(figsize=(5, 5))
plt.errorbar(z_sample, mu_sample, dmu, fmt='.k', ecolor='gray', lw=1)
plt.xlabel(r'$z$')
plt.ylabel(r'$\mu$')
Explanation: Put the data into a Pandas dataframe
End of explanation
snmod_wls = smf.wls("mu ~ z", data=sndata, weights=dmu**-2)
snres_wls = snmod_wls.fit()
print(snres_wls.summary())
z_range = np.arange(0,1.8,0.1)
Explanation: Fit a straight line
End of explanation
fig = plt.figure(figsize=(5, 5))
plt.errorbar(z_sample, mu_sample, dmu, fmt='.k', ecolor='gray', lw=1)
plt.plot(z_sample, snres_wls.fittedvalues, 'b-', label="Weighted LS")
plt.xlabel(r'$z$')
plt.ylabel(r'$\mu$')
Explanation: Overplot the straight line fit
End of explanation
snmod_wls_2nd = smf.wls("mu ~ z + np.power(z,2)", data=sndata, weights=dmu**-2)
Explanation: In the cells below, fit 2nd and 4th degree polynomials and overplot on the data
End of explanation
<END_TASK>