! Metodos Numericos y Simulaciones en Astrofisica
! FCEFN - UNSJ, 2010
!
! ejemplo 4
!
! Aplicacion del metodo bayesiano: sacar bolas de una
! urna con reposicion
!
! Caso de N bolas rojas, M bolas azules, N+M=10
! se sacan en 5*k intentos 3*k rojas y 2*k azules
!
! prob(A|B) es binomial
! prob(B) es uniforme
! prob(A) es una cte. de renormalizacion
!
!
program eje4
implicit none
  integer, parameter ::d=selected_real_kind(16,200)
  integer :: tt,rr,ii
  real(d) :: gammln,comb,prob,rrd,ttd,norm

  write(*,*)'Numero de intentos (>5): '
  read(*,*)tt
  if(tt < 5)stop
  ttd=tt*1._d
  rrd=float(tt*3/5)
  
! calcula el numero combinatorio usando la funcion Gamma 
!
  comb=exp(gammln(ttd+1._d)-gammln(rrd+1._d)-gammln(ttd-rrd+1._d))

! encuentra factor de renormalizacion
!  
  norm=0._d
  do ii=1,10
    norm=norm+comb*(float(ii)/10._d)**rrd*(float(10-ii)/10._d)**(ttd-rrd)
  enddo
  write(*,*)'Const. renormalizacion: ',norm

! calcula prob(B|A) para diferentes valores de N
!
  open(7,file='eje4.out')
  do ii=1,10
    prob=comb*(float(ii)/10._d)**rrd*(float(10-ii)/10._d)**(ttd-rrd)/norm
    write(7,*)ii,prob
  enddo
  close(7)
  
end program eje4

! funcion que devuelve el logaritmo neperiano de la
! funcion gamma.
!
! n! = Gamma(n+1)
!
function gammln(xx)
implicit none
  integer, parameter :: d=selected_real_kind(16,200)
  integer :: jj
  real(d) :: stp,x0,y0,tmp,ser,xx,gammln
  real(d), dimension(6) :: cof

  data cof,stp/76.18009172947146d0,-86.50532032941677d0,24.01409824083091d0, & 
       -1.231739572450155d0,.1208650973866179d-2,-.5395239384953d-5,2.5066282746310005d0/

  x0=xx
  y0=x0
  tmp=x0+5.5_d
  tmp=(x0+0.5_d)*log(tmp)-tmp
  ser=1.000000000190015_d

  do jj=1,6
    y0=y0+1._d
    ser=ser+cof(jj)/y0
  enddo
  gammln=tmp+log(stp*ser/x0)
  return

end
