pro bayesian_rho,r,n, rho,prob_rho,plot=plot,oplot=oplot,pos=pos,$
                 check=check,title0=TITLE0

if(n_params() le 0) then begin
    print,'-----------------------------------------'
    print,'bayesian_rho,r,n,rho,prob_rho,plot=plot,title0=title0'
    print,'-----------------------------------------'
    print,'Calculate the Jeffreys pdf (WJ 4.7) of rho'
    print,'based on data of N pairs with a Pearson r'
    print,'(+assuming a bivariate Gaussians + uniform prior for rho)
    print,'  input: r,n
    print,'  output: rho,prob_rho'
    print,''
    print,'/plot  ->  plot distribution'
    print,'/oplot -> oplot distribution'
    print,' pos   -  returns probability of positive rho'
    print,' check=[nrho,ntrial] -> check the distribution by creating'
    print,'       ntrial bivariate Gaussians for each'
    print,'       nrho values of rho in a range of -1<rho<1'  
    print,'       and sampling the prob of getting r for each rho'
    print,''
    print,'EXAMPLES:'
    print,'bayesian_rho,0.2,50, rho,prob_rho,/plot   ;wj fig4.5 upper right'
    print,'bayesian_rho,0.44,50, rho,prob_rho,/plot   ;wj fig4.5 lower right'
    print,'bayesian_rho,0.2,50, rho,prob_rho,/plot,check=[100,1000] ;slow!'
    print,'-----------------------------------------'
    return
endif

rho=(dindgen(1000)+0.5)/1000.*2.-1.
prob_rho=(1.-rho^2)^(n/2.-0.5)/(1.-rho*r)^(n-1.5)*$
  (1.d0+1.d0/(n-0.5d0)*(1.d0+r*rho)/8.)

;normalize
prob_rho=prob_rho/mean(prob_rho)*.5

;probability of positive rho
pos=mean(prob_rho(where(rho gt 0)))

if(keyword_set(plot)) then begin
    nwin
    title='r='+string(r,'(f6.4)')+'    N='+string(n,'(i6)')
    if(keyword_set(title0)) then title=title0
    plot,rho,prob_rho,xtitle='rho',ytitle='Jeffreys Prob(rho)',$      
      title=title
endif    

if(keyword_set(oplot)) then begin
        oplot,rho,prob_rho
endif    

;---------------------------------------------
;check by creating distribution of prob(rho)
;---------------------------------------------

if(keyword_set(check)) then begin
    
    nrho=check(0)
    ntrial=check(1)

;distribution is checked at these rho values
    rho_num=-1+(findgen(nrho)+0.5)/(nrho-1)*2. ;-1 to 1
    prob_rho_num=rho_num*0.

    for irho=0,nrho-1 do begin

;create ntrial samples for each rho_num
;variances do not affect Pearson r
        rtab=fltarr(ntrial)
        sx=1.
        sy=1.
                
        for i=0l,ntrial-1 do begin
            bivariate_gaussian_f,x,y,n,rho_num(irho),sx=sx,sy=sy
            rtab(i)=correlate(x,y)
        endfor

;estimate the probability of obtaining r calculating the
;fraction of samples giving r at  r \pm 0.02
        ind=where(rtab gt r-0.02 and rtab le r+0.02,count)
        prob_rho_num(irho)=1.*count/ntrial
    endfor

;normalize
    prob_rho_num=prob_rho_num/mean(prob_rho_num)/2.    
    oplot,rho_num,prob_rho_num,psym=10

endif



end