<jupyter_start><jupyter_text>Παρακάτω προσπαθώ να παρουσιάσω την ανάλυση του *Hogg et al 2010 - Data analysis Recipes* η οποία μας λέει από το abstract:Χρησιμοποιούμε πάρα πολλές *considerations involved* ώστε να προσαρμόσουμε ένα μοντέλο σε δεδομένα, όπως για παράδειγμα τη προσαρμογή μιας ευθείας σε σημεία δύο διαστάσεων. Η standard μέθοδος των ελαχίστων τετραγώνων (wls) θεωρεί οτι τα σφάλματα των σημείων στη μια διάσταση είναι *neglible* ενώ στην άλλη περιγράφονται απο Gaussians με γνωστή απόκλιση, υποθέσεις που στη πραγματικότητα συναντώνται εξαιρετικά σπάνια. Θεωρούμε τη περίπτωση ενός ετερογενούς δείγματος δεδομένων με άγνωστα σφάλματα σε δύο διαστάσεις και περιπτώσεις όπου συμπεριλαμβάνονται κακά δεδομένα (outliers) τα οποία εμφανίζουν γραμμική συσχέτιση.*Above all we emphasize the importance of having a “generative model” for the data, even an approx-imate one. Once there is a generative model, the subsequent fittingis non-arbitrary because the model permits direct computation of thelikelihood of the parameters or the posterior probability distribution.Construction of a posterior probability distribution is indispensible ifthere are “nuisance parameters” to marginalize away.*<jupyter_code>import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import rc,rcParams
rc('text', usetex=True)
rcParams['figure.figsize'] = (8, 6.5)
rcParams['ytick.labelsize'],rcParams['xtick.labelsize'] = 17.,17.
rcParams['axes.labelsize']=19.
rcParams['legend.fontsize']=17.
rcParams['axes.titlesize']=20.
rcParams['text.latex.preamble'] = ['\\usepackage{siunitx}']
import seaborn
seaborn.despine()
seaborn.set_style('white', {'xes.linewidth': 0.5, 'axes.edgecolor':'black'})
seaborn.despine(left=True)
%load_ext autoreload
from astroML.linear_model import TLS_logL
from astroML.plotting.mcmc import convert_to_stdev
from astroML.datasets import fetch_hogg2010test
from astropy.table import Table
from scipy import polyfit,linalg
from scipy.optimize import curve_fit
import scipy.stats as sts<jupyter_output>The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload<jupyter_text>Η βιβλιοθήκη astroML (που μπορει να εγκατασταθεί έυκολα μέσω της ```sudo pip install astroml``` εμπεριέχει τα δεδομένα που χρησιμοποιούνται στο paper. Παρόλαυτα όποιος τα θέλει αλλιώς ή δεν θέλει/μπορεί να εγκαταστήσει τη βιβλιοθήκη μπορεί να μου τα ζητήσει.<jupyter_code>data = fetch_hogg2010test()
data = data[5:]  # no outliers
x = data['x']
y = data['y']
sigma_x = data['sigma_x']
sigma_y = data['sigma_y']
rho_xy = data['rho_xy']
np.savetxt('data',data)
t = Table([x, y, sigma_x,sigma_y,rho_xy], names=('X', 'Y', '$\sigma _x$','$\sigma _y$','$\\rho _{xy}$'))#, meta={'name': 'first table'})
t
plt.errorbar(x,y,yerr=sigma_y,xerr=sigma_x,fmt='o')<jupyter_output><empty_output><jupyter_text>Έχουμε ένα σύνολο $N>2$ σημείων $(x_i,y_i)$ με γνωστές gaussian αποκλίσεις $\sigma _{yi}$ στη $y$ διεύθυνση και χωρίς καμία απροσδιοριστία στη $x$ διεύθυνση.Αναζητούμε τη συνάρτηση $f(x)$ της μορφής: $ f(x)=\theta _1 x+\theta _0 $ η οποία προσαρμόζεται καλύτερα στα σημεία.Κατασκευάζουμε τους πίνακες \begin{align}\mathbf{Y}&=\begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N\end{pmatrix} \\\mathbf{A}&=\begin{pmatrix} 1 & x_1 \\ 1 & x_2 \\ \vdots & \vdots \\ 1 & x_N \end{pmatrix} \\\mathbf{C}&=\begin{pmatrix} \frac{1}{\sigma ^2_{y1}} & 0 & 0 & \cdots & 0 \\ 0 & \frac{1}{\sigma ^2_{y2}} & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & \frac{1}{\sigma ^2_{yN}} \end{pmatrix}\end{align}Η λύση καλύτερης προσαρμογης δίνεται από τη:\begin{equation}\begin{pmatrix} b \\ m \end{pmatrix} = \mathbf{S}= [\mathbf{A}^T \mathbf{C} \mathbf{A}]^{-1} [\mathbf{A}^T \mathbf{C} \mathbf{Y}]\end{equation}<jupyter_code>Y=y
A=np.stack((np.ones(x.shape[0]),x),axis=-1)
C=np.diagflat(1./sigma_y**2)
first=np.linalg.inv(A.T.dot(C.dot(A)))
second=A.T.dot(C.dot(Y))
M=np.dot(first,second)
print(M)<jupyter_output>[34.84594131  2.23253622]<jupyter_text>***Ερώτηση:*** Πως μπορούμε να υπολογίσουμε τα σφάλματα μέσω αναλυτικής έκφρασης?<jupyter_code>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)))
b_ls,m_ls<jupyter_output><empty_output><jupyter_text>Εκτίμηση των παραμέτρων προσαρμογής μέσω optimizationΜέσω της βιβλιοθήκης ```scipy.optimization```  Εκτίμηση των παραμέτρων προσαρμογής μέσω της μεθόδου ελαχίστων τετραγώνωνΜέσω της βιβλιοθήκης ```scipy.curve_fit```<jupyter_code>def f(x,b,m): return b+m*x
Mcf,dMcf2=curve_fit(f,x,y,[2.,30.],sigma=sigma_y)
dMcf=np.sqrt(np.diag(dMcf2))
print(Mcf, dMcf)
print('Y=({:.1f} +/- {:.1f}) + ({:.1f} +/- {:.2f}) X'.format(Mcf[0],dMcf[0],Mcf[1],dMcf[1]))
#----PLOT----------------
xx=np.linspace(x.min(),x.max())
plt.errorbar(x,y,yerr=sigma_y,fmt='o')
plt.plot(xx,f(xx,Mcf[0],Mcf[1]),label='Scipy Curve Fit')
plt.fill_between(xx,f(xx,Mcf[0]-2.*dMcf[0],Mcf[1]-2.*dMcf[1]),f(xx,Mcf[0]+2.*dMcf[0],Mcf[1]+2.*dMcf[1]),alpha=0.25,label='95\% Confidence')
#plt.fill_between(xx,f(xx,Mcf[0]-2.*dMcf[0],Mcf[1]),f(xx,Mcf[0]+2.*dMcf[0],Mcf[1]),alpha=0.25,label='95\% Confidence')
#plt.fill_between(xx,f(xx,Mcf[0],Mcf[1]-2.*dMcf[1]),f(xx,Mcf[0],Mcf[1]+2.*dMcf[1]),alpha=0.25,label='95\% Confidence')
plt.legend(loc='best')
N=len(x)
d=N-2
NN=200
bb=np.linspace(-40,110,NN)
mm=np.linspace(1.75,2.7,NN)
LogL=np.zeros((len(bb),len(mm)))
for i,b in enumerate(bb):
    for j,m in enumerate(mm):
        LogL[i,j]=np.sum((y-f(x,b,m))**2/sigma_y**2)
print(LogL.min(),LogL.min()/d,(LogL.min()/d-1.)/np.sqrt(2./d))
ax=plt.subplot()
c=ax.contour(LogL,np.linspace(LogL.min(),1.8*LogL.min(),5),extent=[mm.min(),mm.max(),bb.min(),bb.max()])
c2=ax.contourf(LogL,[LogL.min(),LogL.min()+1.],extent=[mm.min(),mm.max(),bb.min(),bb.max()],label='$\chi ^2+1$')
ax.set_xlabel('slope')
ax.set_ylabel('intercept')
ax.vlines(Mcf[1],bb.min(),bb.max(),linestyles='--',alpha=0.2)
ax.hlines(Mcf[0],mm.min(),mm.max(),linestyles='--',alpha=0.2)
ax.vlines(Mcf[1]+dMcf[1],Mcf[0],Mcf[0]+dMcf[0],linestyles='--',alpha=0.7)#,label=r'Intercept $\sigma$')
ax.hlines(Mcf[0]-dMcf[0],Mcf[1],Mcf[1]-dMcf[1],linestyles='--',alpha=0.7)#,label=r'Slope $\sigma$')
ax.plot([Mcf[1]],[Mcf[0]],'o',label='Minimum at [{:.2f},{:.2f}]'.format(Mcf[1],Mcf[0]))
ax.set_title('$\chi ^2$ Contours')
plt.colorbar(c,ax=ax)
ax.legend()<jupyter_output>18.585270532879612 1.4296361948368932 1.095361670609535<jupyter_text>Εδώ προσπαθώ να υπολογίσω τη παράμετρο p. Η παράμετρος p ουσιαστικά είναι η πιθανότητα η κάθε παράμετρος να μην παίζει ρόλο. Δηλαδή αν για παράδειγμα η ευθεία ειχε κλίση αλλα το σφάλμα της κλίση είναι μέγαλο τότε υπάρχει μια πιθανότητα η κλιση να είναι 0 άρα ουσιαστικά τα δεδομένα να μην έχουν αυτο το βαθμό ελευθερίας που εμεις θεωρούμε.Τέλος Πάντων είναι λίγο άσχετο με την όλη θεματολογια και μπορεί να το βγάλω στο τέλος.<jupyter_code>from scipy.stats import chisquare
y,f(x,Mcf[0],Mcf[1])
chisquare(f_obs=y,f_exp=f(x,Mcf[0],Mcf[1]),ddof=1)
s=np.random.standard_t(d, size=100000)
tt=np.array([Mcf[0]/dMcf[0]/np.sqrt(N),Mcf[1]/dMcf[1]/np.sqrt(N)])
print np.sum(s>tt[0])/float(len(s))+np.sum(s<-tt[0])/float(len(s))
print np.sum(s>tt[1])/float(len(s))+np.sum(s<-tt[1])/float(len(s))
def f2(x,b,m1,m2): return b+m1*x+m2*x**2
M2cf,dM2cf2=curve_fit(f2,x,y,[2.,30.,2],sigma=sigma_y)
dM2cf=np.sqrt(np.diag(dM2cf2))
print(M2cf, dM2cf)
print('Y=({:.1f} +- {:.1f}) + ({:.1f} +- {:.2f}) X+({:.3f} +- {:.4f}) X^2'.format(M2cf[0],dM2cf[0],M2cf[1],dM2cf[1],M2cf[2],dM2cf[2]))
xx=np.linspace(x.min(),x.max())
plt.errorbar(x,y,yerr=sigma_y,fmt='o')
plt.plot(xx,f2(xx,M2cf[0],M2cf[1],M2cf[2]),label='Scipy Curve Fit')
plt.fill_between(xx,f2(xx,M2cf[0]-2.*dM2cf[0],M2cf[1]-2.*dM2cf[1],M2cf[2]-2.*dM2cf[2]),
                 f2(xx,M2cf[0]+2.*dM2cf[0],M2cf[1]+2.*dM2cf[1],M2cf[2]+2.*dM2cf[2]),
                 alpha=0.25,label='95\% Confidence')
plt.legend(loc='best')
#plt.ylim([0,700])
tt2
N=x.shape[0]
d=N-3
s=np.random.standard_t(d, size=100000)
tt2=np.array([M2cf[0]/dM2cf[0]/np.sqrt(N),M2cf[1]/dM2cf[1]/np.sqrt(N),M2cf[2]/dM2cf[2]/np.sqrt(N)])
plt.hist(s, bins=100, normed=True)
plt.vlines(tt2[2],0,0.4)
print(tt2)
print(np.sum(s>tt2[0])/float(len(s))+np.sum(s<-tt2[0])/float(len(s)))
print(np.sum(s>tt2[1])/float(len(s))+np.sum(s<-tt2[1])/float(len(s)))
print(np.sum(s>tt2[2])/float(len(s))+np.sum(s<-tt2[2])/float(len(s)))<jupyter_output>[0.39977454 0.59171963 0.23587458]
0.69722
0.56693
0.81594<jupyter_text>Data with outliers<jupyter_code>data = fetch_hogg2010test()
x = data['x']
y = data['y']
sigma_x = data['sigma_x']
sigma_y = data['sigma_y']
rho_xy = data['rho_xy']
plt.errorbar(x,y,yerr=sigma_y,xerr=sigma_x,fmt='o')
Mcft,dMcf2t=curve_fit(f,x,y,[2.,30.],sigma=sigma_y)
dMcft=np.sqrt(np.diag(dMcf2t))
print Mcft, dMcft
print 'Y=({:.1f} +- {:.1f}) + ({:.1f} +- {:.2f}) X'.format(Mcft[0],dMcft[0],Mcft[1],dMcft[1])
xx=np.linspace(x.min(),x.max())
plt.errorbar(x,y,yerr=sigma_y,fmt='o')
plt.plot(xx,f(xx,Mcft[0],Mcft[1]),label='Scipy Curve Fit')
plt.fill_between(xx,f(xx,Mcft[0]-2.*dMcft[0],Mcft[1]),f(xx,Mcft[0]+2.*dMcft[0],Mcft[1]),color='g',alpha=0.3,label='95\% Confidence')
plt.fill_between(xx,f(xx,Mcft[0],Mcft[1]-2.*dMcft[1]),f(xx,Mcft[0],Mcft[1]+2.*dMcft[1]),color='g',alpha=0.3,label='95\% Confidence')
plt.plot(xx,f(xx,Mcf[0],Mcf[1]),'--',label='Scipy Curve Fit (No outliers)')
plt.legend(loc='best')
N=200
bb=np.linspace(-40,320,N)
mm=np.linspace(0.5,2.4,N)
LogL=np.zeros((len(bb),len(mm)))
for i,b in enumerate(bb):
    for j,m in enumerate(mm):
        LogL[i,j]=np.sum((y-f(x,b,m))**2/sigma_y**2)
print LogL.min(),LogL.mean()
ax=plt.subplot()
c=ax.contour(LogL,np.linspace(LogL.min(),1.8*LogL.min(),5),extent=[mm.min(),mm.max(),bb.min(),bb.max()])
c2=ax.contourf(LogL,[LogL.min(),LogL.min()+5],extent=[mm.min(),mm.max(),bb.min(),bb.max()],label='$\chi ^2+1$')
ax.set_xlabel('slope')
ax.set_ylabel('intercept')
ax.vlines(Mcft[1],bb.min(),bb.max(),linestyles='--',alpha=0.2)
ax.hlines(Mcft[0],mm.min(),mm.max(),linestyles='--',alpha=0.2)
ax.vlines(Mcft[1]+dMcft[1],Mcft[0],Mcft[0]+dMcft[0],linestyles='--',alpha=0.7)#,label=r'Intercept $\sigma$')
ax.hlines(Mcft[0]-dMcft[0],Mcft[1],Mcft[1]-dMcft[1],linestyles='--',alpha=0.7)#,label=r'Slope $\sigma$')
ax.plot([Mcft[1]],[Mcft[0]],'o',label='Minimum at [{:.2f},{:.2f}]'.format(Mcft[1],Mcft[0]))
ax.set_title('$\chi ^2$ Contours')
plt.colorbar(c,ax=ax)
ax.legend()
def f2(x,b,m1,m2): return b+m1*x+m2*x**2
M2cf,dM2cf2=curve_fit(f2,x,y,[2.,30.,2],sigma=sigma_y)
dM2cf=np.sqrt(np.diag(dM2cf2))
print M2cf, dM2cf
print 'Y=({:.1f} +- {:.1f}) + ({:.1f} +- {:.2f}) X+({:.3f} +- {:.4f}) X^2'.format(M2cf[0],dM2cf[0],M2cf[1],dM2cf[1],M2cf[2],dM2cf[2])
xx=np.linspace(x.min(),x.max())
plt.errorbar(x,y,yerr=sigma_y,fmt='o')
plt.plot(xx,f2(xx,M2cf[0],M2cf[1],M2cf[2]),label='Scipy Curve Fit')
plt.fill_between(xx,f2(xx,M2cf[0]-2.*dM2cf[0],M2cf[1],M2cf[2]),f2(xx,M2cf[0]+2.*dM2cf[0],M2cf[1],M2cf[2]),alpha=0.25,label='95\% Confidence')
plt.fill_between(xx,f2(xx,M2cf[0],M2cf[1]-2.*dM2cf[1],M2cf[2]),f2(xx,M2cf[0],M2cf[1]+2.*dM2cf[1],M2cf[2]),alpha=0.25,label='95\% Confidence')
plt.fill_between(xx,f2(xx,M2cf[0],M2cf[1],M2cf[2]-2.*dM2cf[2]),f2(xx,M2cf[0],M2cf[1],M2cf[2]+2.*dM2cf[2]),alpha=0.25,label='95\% Confidence')
plt.legend(loc='best')
#plt.ylim([0,700])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.22985233944138053, 'fraction_numerical': 0.057818893435331795, 'mean_word_length': 7.111832611832612, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '226350', 'license': 'MIT', 'n_tokens_mistral': 6708, 'n_tokens_neox': 5695, 'n_words': 739, 'path': 'Linear Regression.ipynb', 'repo_name': 'Mixpap/astrostatistics', 'size': '11242'}