; scripts to do exercise 3 // phy 535 // fall 2011
;
; instructions
;	at the IDL prompt type
;
;	IDL> .run soln_exercise3
;

window,0
; read in data
readcol,'exercise3.dat',x,sigmax,y,sigmay,z,sigmaz
print,' '

; do l-s fit -- here are the equations
delta = n_elements(x) * total(x*x) - total(x)^2
bb = ( ( total(x*x) * total(y) ) - ( total(x) * total(x*y) ) ) / delta
mm = ( n_elements(x) * total(x*y) - total(x) * total(y) ) / delta
print,'unweighted L/S fit, from equations: ',bb,mm

; now i compare to a canned routine to see how i did
result = poly_fit(x,y,1,chisq=csunw)

; result[0] is the y intercept (y value when x=0), result[1] is the slope
print,'unweighted L/S fit, from canned routine: ',result[0],result[1]

; plot data and best fit line
ploterror,x,y,sigmax,sigmay,psym=4,xrange=[-3,20],yrange=[-10,120],xstyle=1,ystyle=1
oplot,findgen(12)-1.,bb+mm*(findgen(12)-1.)

; now do weighted fit - here are the equations
; let 
; 1/w=sigma_x^2 for the x only terms, 
; 1/w=sigma_y^2 for the y only terms, and
; 1/w=sigma_x*sigma_y for the cross terms
; i'm also going to let 1/w=sigma_x*sigma_y for the no-x, no-y terms, though
;	there are other choices that are possible
wx = 1./(sigmax*sigmax)
wy = 1./(sigmay*sigmay)
wxy = 1./(sigmax*sigmay)
w = 1./(sigmax*sigmay)
deltaw = total(w) * total(wx*x*x) - total(wx*x)^2
bbw = ( ( total(wx*x*x) * total(wy*y) ) - ( total(wx*x) * total(wxy*x*y) ) ) / deltaw
mmw = ( total(w) * total(wxy*x*y) - total(wx*x) * total(wy*y) ) / deltaw
print,'weighted L/S fit, from equations: ',bbw,mmw
; plot this line
oplot,findgen(12)-1.,bbw+mmw*(findgen(12)-1.),linestyle=1

; now i compare to a canned routine to see how i did
; note that this canned routine can only accept errors
;	on y -- another reason to not use canned routines
;	but use the full on equations instead
result = poly_fit(x,y,1,measure_errors=sigmay,chisq=csw)
print,'weighted L/S fit, from canned routine: ',result[0],result[1]
; plot this line
oplot,findgen(12)-1.,result[0]+result[1]*(findgen(12)-1.),linestyle=2

print,' '

; print chisq value
print,'unweighted chisq: ',csunw/(10.-2.)
; reduced chisq is ~9.5

print,'weighted chisq: ',csw/(10.-2.)
; reduced chisq is ~0.03

; now we need to calculate covariance and correlation
meanx = mean(x)
momentx = moment(x)
meany = mean(y)
momenty = moment(y)

print,' '
; now calculate covariance
print,'covariance: ',(total((x-meanx)*(y-meany)))/10.
; value is around 41

; let's check. sqrt(momentx[1]) is sigmax, sqrt(momenty[1]) is sigmay
print,'sanity check:'
print,'---> [sigma_x, sigma_y, sigma_x * sigma_y]: ',sqrt(momentx[1]),sqrt(momenty[1]),sqrt(momentx[1])*sqrt(momenty[1])
print,'---> covariance is a little less than sigma_x * sigma_y -- so strongly covary, but not as unity'
print,' ' 

; correlation:
ssxy = total( (x - meanx)*(y - meany) ) 
ssxx = total( (x - meanx)*(x - meanx) )
ssyy = total( (y - meany)*(y - meany) )
rr = ssxy/( sqrt(ssxx) * sqrt(ssyy) )
print,'correlation coefficient: ',rr
; i get 0.98

; what is the probability that x and y are correlated? see lookup table
print,'---> see the lookup table for the probability that x and y are correlated'
print,' '

print,' '
print,'NOTE: type .continue to proceed to question 2'
print,' '
stop

;;;;;;;;
; now let's do question 2
;
;

window,0
loadct,13
ploterror,x,y,sigmax,sigmay,psym=4,xrange=[-3,20],yrange=[-10,120],xstyle=1,ystyle=1
oploterror,x,z,sigmax,sigmaz,psym=4,color=200,errcolor=200
loadct,0
window,1

; then do t test and K-S test (see exercise2.pro for results here)
print,' '

; t test 
; mean of y is meany
; moment of y is momenty
meanz = mean(z)
momentz = moment(z)
df = 18. ; Na - 1 + Nb - 1 = 18
; sigma_y^2 = variance(y) = momenty[1]
sigmap = sqrt(((10.*momenty[1]) + (10.*momentz[1])) / df)
sigmamm = sqrt(((sigmap^2)/10.)+((sigmap^2)/10.))
tt = (meany - meanz) / sigmamm
print,'t = ',tt,' ; df = ',df
print,'---> see the lookup table'
print,' '

; K-S test
sorty = y(sort(y))
sortz = z(sort(z))
plot,sorty,(findgen(10)+1.)/10.,psym=-4,title='K-S test',xrange=[-10,60],xstyle=1
oplot,sortz,(findgen(10)+1.)/10.,psym=-5

print,' '
print,'K-S test: find distance D and use the lookup table'
print,' '
print,'check against canned routine:'

; check against canned routine
kstwo,y,z,D,prob

print,'KS statistic: ',D
print,'prob two data sets are the same: ',prob

print,' '
end