;---------------------------------------------------
;Illustrate the Bayesian (Jeffrey's) pdf of rho
;based on having data of N pairs with Pearson r.
; HS 041105
;---------------------------------------------------


program='bayesian_rho_driver'
ps=0.1

;goto,plotb


print,'TAMA KESTAA TOVIN ...'
;--------------------------------------------------------------
;r varied
;--------------------------------------------------------------
psdirect,program,ps,color=1

!p.linestyle=2
bayesian_rho,0.75,10,rho,prob_rho,/plot,$
title='Different values of r, for samples with N=10'
!p.linestyle=0 & xyouts,0.75,max(prob_rho)+0.05,ali=0.5,'r=0.75'

thick,3
!p.linestyle=0
bayesian_rho,0.5,10,rho,prob_rho,/oplot,check=[50,5000]  ;THIS TAKES AWHILE!
!p.linestyle=0 & xyouts,0.5,max(prob_rho)+0.1,ali=0.5,'r=0.50'
thick,1

!p.linestyle=1
bayesian_rho,0.25,10,rho,prob_rho,/oplot
!p.linestyle=0 & xyouts,0.25,max(prob_rho)+0.05,ali=0.5,'r=0.25'

!p.linestyle=1
bayesian_rho,-0.25,10,rho,prob_rho,/oplot
!p.linestyle=0 & xyouts,-0.25,max(prob_rho)+0.05,ali=0.5,'r=-0.25'

!p.linestyle=0
bayesian_rho,-0.50,10,rho,prob_rho,/oplot
!p.linestyle=0 & xyouts,-0.5,max(prob_rho)+0.05,ali=0.5,'r=-0.50'

!p.linestyle=2
bayesian_rho,-0.75,10,rho,prob_rho,/oplot
!p.linestyle=0 & xyouts,-0.75,max(prob_rho)+0.05,ali=0.5,'r=-0.75'

;-------

bayesian_rho,0,10,rho,prob_rho,/oplot
!p.linestyle=0 & xyouts,0,max(prob_rho)+0.05,ali=0.5,'r=0'
xyouts,-0.1,2.2,'r=0.5 verified numerically !c!c  with Nrho=50, Ntrial=5000'

 
defplot
psdirect,program,ps,color=1,/stop

;--------------------------------------------------------------
;N varied
;--------------------------------------------------------------

plotb:   
psdirect,program+'_b',ps,color=1


!p.linestyle=0
bayesian_rho,0.25,1000,rho,prob_rho,/plot,$
title='Fixed r=0.25, samples with different N'
!p.linestyle=0 & xyouts,0.5,max(prob_rho)+0.05,ali=0.5,'N=1000'

!p.linestyle=2
bayesian_rho,0.25,100,rho,prob_rho,/oplot
!p.linestyle=0 & xyouts,0.5,max(prob_rho)+0.05,ali=0.5,'N=100'

!p.linestyle=1
bayesian_rho,0.25,10,rho,prob_rho,/oplot
!p.linestyle=0 & xyouts,0.5,max(prob_rho)+0.05,ali=0.5,'N=10'

psdirect,program+'_b',ps,/stop

nwin

;--------------------------------------------------------------
;Probability of positive rho vs r and N
;--------------------------------------------------------------

plotc:   
psdirect,program+'_c',ps,color=1,xsize=12,ysize=7


n100=200
ntab=lindgen(n100)+1
postab25=findgen(n100)
postab50=findgen(n100)
postab75=findgen(n100)
postab15=findgen(n100)
postab10=findgen(n100)

for i=0,n_elements(ntab)-1 do begin
    bayesian_rho,0.10,ntab(i),rho,prob_rho,pos=pos
    postab10(i)=pos
    bayesian_rho,0.15,ntab(i),rho,prob_rho,pos=pos
    postab15(i)=pos
    bayesian_rho,0.25,ntab(i),rho,prob_rho,pos=pos
    postab25(i)=pos
    bayesian_rho,0.50,ntab(i),rho,prob_rho,pos=pos
    postab50(i)=pos
    bayesian_rho,0.75,ntab(i),rho,prob_rho,pos=pos
    postab75(i)=pos
endfor


!p.multi=[0,1,1]

plot,ntab,postab75,xr=[10,n100],xs=1,yr=[0.89,1.01],ys=1,$
  xtitle='N',ytitle='Prob rho>0'
oplot,ntab,postab50,lines=2
oplot,ntab,postab25,lines=1
oplot,ntab,postab15,lines=3
oplot,ntab,postab10,lines=4

xyouts,150,0.90,'r=0.10'
xyouts,85,0.90,'r=0.15'
xyouts,40,0.92,'r=0.25'
xyouts,17,0.96,'r=0.50'
xyouts,12,0.985,'r=0.75'

oplot,[10,n100],[1,1]*.95,lines=1

psdirect,program+'_c',ps,/stop

end