;REDO THE LEAST SQUARES FIT, DOING A CHI SQUARES FIT INSTEAD... ;follow notation in section 8.1.4... ;SUBSCRIPT ALL CHI FIT ARRAYS BY _CHI... x = dblarr(3,4) m=4 n=3 ;FROM THE LEAST SQUARES FIT, WE FIND... sigma = sqrt( 20.d0) time = double([5,7,9,11]) ;time=time-8. x[0,*]=1. x[1,*]=time x[2,*]=time^2 y= double([142,168,211,251]) ydata = y y = transpose( y) ;DEFINE THE SIGMAS... sigma = fltarr(M) + 1. sigma = dblarr(M) + sqrt(20.d0) ;sigma = 1. + findgen(m) ;CREATE A DIAGONAL MATRIX WHOSE ELEMENTS ARE THE INVERSES OF THE SIGMAS... w = fltarr( m,m) w [ indgen(m) * (m+1)] = 1./sigma ;DERIVE WEIGHTED VERSIONS OF X AND Y xw = w ## x yw = w ## y xxw= transpose( xw) ## xw xyw = transpose(xw) ## yw xxwi = invert(xxw) a = xxwi ## xyw ybarw = xw ## a delyw = yw - ybarw ybar = invert( w) ## xw ## a dely = y - ybar sigmasq = total( dely^2)/(m-n) chisq = total((yw - ybarw)^2) reduced_chisq = total((yw - ybarw)^2)/(m-n) dcw = xxwi[ (n+1) * indgen(n)] vardc_official = dcw vardc_realworld = reduced_chisq* dcw sigsqdc_official = sqrt( vardc_official) sigsqdc_realworld = sqrt( vardc_realworld) ncov_chi = xxwi/sqrt(dcw # dcw) ;--------------------------------------------------- ;stop ;---------------- indx_new = [0,1] ;indx_new = [0] indx_new= [1,2] xxwi_new = pminor( indx_new, xxwi) xxw_new = invert( xxwi_new) ;--------------- do the case of getting uncerts in a0 alone ------ ;PICK VALUES FOR A1 AND A2. SOLVE FOR A0. DO WE GET THE CORRECT ;RESULT FOR A0, WHICH IS THE WEIGHTED AVG OF YDAT? ;answer: for the average a_n, yes. ; the sample variance is still 20. the variance of the mean is 20./4. ; so the error of the mean is sqrt(20/4)=2.236. so xxwi_n should equal ; 5. AND IT IS!! ;for the variance, which is xxwi_n, we should get the sigma^2 of ;the points. ;now do ls fit for a_0 using the derived values for ydat = ydata - a[1]*time - a[2]*time^2 ;ydat = ydata - (a[1]+9.)*time - (a[2]-.559)*time^2 x_n = dblarr( 1, 4) x_n[ 0,*] = 1. xw_n = w ## x_n yw_n = w ## ydat xxw_n = transpose( xw_n) ## xw_n xyw_n = transpose( xw_n) ## yw_n xxwi_n = invert( xxw_n) a_n = xxwi_n ## xyw_n ybarw_n = xw_n ## a_n delyw_n = yw_n - ybarw_n ybar_n = invert( w) ## xw_n ## a_n dely_n = ydat - ybar_n sigmasq = total( dely_n^2)/(m-n) chisq_n = total( delyw_n^2) vardc_n = xxwi_n[ 0] sigdc_n = sqrt( vardc_n) ;now the question is: if we perturb (a[1],1[2]) by the proper amount, ;does delta chisq equal 1? go back and redefine ydat, substituting ;a[1] + 9., a[2] = .559. THIS WORKS.. ;conclusion: we derived delta chisq equals 1 in two ways. one using the ;curvature matrix below. the other doing a separate chi squared fit above. ;thus the two methods are equivalent. ;------------------------------------- stop polyfit, time, ydat, 0, coeffs, sigcoeffs, yfit, sigma stop dela0_v = 34./100.*(findgen( 400) - 200.) dela0_v = 9./100.*(findgen( 400) - 200.) dela1_v = 0.6/100.*(findgen( 400) - 200.) dela01 = fltarr( 400,400) for nr=0,399 do dela01[*, nr]= dela0_v for nr=0,399 do dela01[nr, *]= dela1_v delchisq = fltarr( 400,400) ;stop for nrx=0,399 do begin &$ for nry=0, 399 do begin &$ dela = [dela0_v[ nrx], dela1_v[ nry]] &$ delchisq[ nrx, nry] = $ dela ## xxw_new ## transpose(dela) &$ endfor &$ endfor levels = [1,2.3] contour, delchisq, dela0_v, dela1_v, levels=levels plots, /data, 34,-9,psym=4 plots, sigsqdc_official[ indx_new[0]], sigsqdc_official[ indx_new[1]],psym=2 plots, -sigsqdc_official[ indx_new[0]], sigsqdc_official[ indx_new[1]],psym=2 plots, sigsqdc_official[ indx_new[0]], -sigsqdc_official[ indx_new[1]],psym=2