;----------------------------------------------------------
;illustrate the error estimates of linear least-squares fit
;----------------------------------------------------------
program0='linear_fit_example'

  ps=0.1

  for icase=1,3 do begin


;normal case
      if(icase eq 1) then begin
          addtitle='GAUSSIAN ERRORS'
          program=program0+'_a'
      endif

;check the effect of assuming gaussian errors
;use uniform errors with the same variance (sqrt(12)...)
      if(icase eq 2) then begin
          addtitle='UNIFORM ERRORS'
          program=program0+'_b'
      endif

      if(icase eq 3) then begin
          addtitle='EXPONENTIAL ERRORS'
          program=program0+'_c'
      endif

  psdirect,program,ps
;create data according to a linear dependence Y=ax+b
;+ Gaussian error 

  a0=1
  b0=1

  nx=50
  x=findgen(nx)
  std=x*0.+5*2

  seed=2
  if(icase eq 1) then y=a0*x+b0+randomn(seed,nx)*std
  if(icase eq 2) then y=a0*x+b0+(randomu(seed,nx)-0.5)*std*sqrt(12)
  if(icase eq 3) then y=a0*x+b0+randomu(seed,nx,gamma=1)*std*sign(randomu(seed,nx)-0.5)/sqrt(2.)

;note that in Numerical Recipes standard form is y=a+bx
;but we stick with Wall-Jenkins notation y=ax+b
  ba=linfit(x,y,measure_errors=std,cov=cov,chisq=chisq,prob=prob)
  b=ba(0)
  a=ba(1)
;cov= 2*2 covariance matrix estimated as explained in the lectures


  !P.MULTI=[0,2,2]
  nwin
  plot,x,y,psym=6,col=2
  oplot,x,a*x+b,col=2,thick=3

  print,'----------------------------------------------------'
  print,'linear_fit_example.pro
  print,'----------------------------------------------------'
  print,'true:   a,b',a0,b0
  print,'fitted: a,b',a,b
  print,'chi2',chisq
  print,'prob (>chi2)',PROB
  print,'sigma_a',sqrt(cov(1,1)) ;note the meaning of a,b!
  print,'sigma_b',sqrt(cov(0,0)) ;note the meaning of a,b!
  print,'cov(a,b)',cov(0,1)

  
;compare these error estimates with that obtained
;by making 10000 random realizations of the x,y relation


  plot,[0,1],[0,1],/nod,xs=15,ys=15   ; empty
  plot,x,y,psym=6,col=2
  oplot,x,a*x+b,col=2,thick=3

  atab=fltarr(1000)
  btab=atab
  seed=1
  for i=0,1000-1 do begin
      if(icase eq 1) then  y=a0*x+b0+randomn(seed,nx)*std
      if(icase eq 2) then  y=a0*x+b0+(randomu(seed,nx)-0.5)*std*sqrt(12)
      if(icase eq 3) then  y=a0*x+b0+randomu(seed,nx,gamma=1)*std*sign(randomu(seed,nx)-0.5)/sqrt(2.)
      ba=linfit(x,y,measure_errors=std)
      atab(i)=ba(1)
      btab(i)=ba(0)
if i mod 20 eq 0 then begin
      oplot,x,x*ba(1)+ba(0),lines=2
      oplot,x,y,psym=3
      endif
  endfor


  plot,atab,btab,psym=3,XTITLE='a',ytitle='b',xr=[0.6,1.4],xs=1,yr=[-15,15],ys=1
  plots,a0,b0,psym=1,col=2,syms=3

  var_a=variance(atab)
  var_b=variance(btab)
  cov_ab=correlate(atab,btab,/cov)

  print,'----------------------------------------------------'
  print,'1000 TRIALS:'
  print,'sigma_a',sqrt(VAR_A)  
  print,'sigma_b',sqrt(VAR_B)  
  print,'cov(a,b)',COV_AB
  print,'----------------------------------------------------'

  xp=0.53
  yp=0.95
  dy=-0.025

  xyouts,/nor,xp,yp,'-----------------------------------'& yp=yp+dy
  xyouts,/nor,xp,yp,'linear_fit_example.pro '+addtitle   & yp=yp+dy 
  xyouts,/nor,xp,yp,'-----------------------------------'& yp=yp+dy
  xyouts,/nor,xp,yp,'true:   a,b'+string(a0)+string(b0)  & yp=yp+dy
  xyouts,/nor,xp,yp,'fitted: a,b'+string(a)+string(b)    & yp=yp+dy
  xyouts,/nor,xp,yp,'chi2'+string(chisq)                 & yp=yp+dy
  xyouts,/nor,xp,yp,'prob (>chi2)'+string(PROB)          & yp=yp+dy
  xyouts,/nor,xp,yp,'sigma_a'+string(sqrt(cov(1,1)))     & yp=yp+dy
  xyouts,/nor,xp,yp,'sigma_b'+string(sqrt(cov(0,0)))     & yp=yp+dy
  xyouts,/nor,xp,yp,'cov(a,b)'+string(cov(0,1))          & yp=yp+dy
  xyouts,/nor,xp,yp,'----------------------------------' & yp=yp+dy
  xyouts,/nor,xp,yp,'1000  TRIALS:'                      & yp=yp+dy
  xyouts,/nor,xp,yp,'sigma_a'+string(sqrt(VAR_A))        & yp=yp+dy
  xyouts,/nor,xp,yp,'sigma_b'+string(sqrt(VAR_B))        & yp=yp+dy
  xyouts,/nor,xp,yp,'cov(a,b)'+string(COV_AB)            & yp=yp+dy
  xyouts,/nor,xp,yp,'-----------------------------------'& yp=yp+dy
 

  !P.MULTI=0
  psdirect,program,ps,/stop
  endfor


;check the errors used

  psdirect,program+'_d',ps


  nx=100000l
  std=1.
  for icase=1,3 do begin
      if(icase eq 1) then  y=randomn(seed,nx)*std
      if(icase eq 2) then  y=(randomu(seed,nx)-0.5)*std*sqrt(12)
      if(icase eq 3) then  y=randomu(seed,nx,gamma=1)*std*sign(randomu(seed,nx)-0.5)/sqrt(2.)
      
      histo_f,y,-5,5,.1,xx,yy,gg
      
      if(icase eq 1) then plot,xx,yy,/ylog,yr=[0.0001,1],xtitle='x/std',ytitle='log(pdf)'
      if(icase ne 1) then oplot,xx,yy,col=icase,lines=icase
  endfor
  
label_data,0.35,0.3,['gaussian','uniform','exponential'],lines=[0,2,3]
psdirect,program+'_d',ps,/stop

end