! Metodos Numericos y Simulaciones en Astrofisica
! FCEFN - UNSJ, 2010
!
! ejemplo 4a
!
! Aplicacion del metodo bayesiano: fuente de un cierto flujo
! observado
! 
! El flujo observado es f y corresponde a un flujo verdadero
! F con distribucion normal y variancia s^2=1. La probabilidad 
! a priori es arbitraria y F se encuentra en el rango [1:4]
!
! Los valores observados para f son:
!
! [2., 1.3, 3., 1.5, 2., 1.8]
!
! prob(A|B) es normal
! prob(B) es arbitraria: prob(B)=K_1 * F^(P)
! prob(A) es una cte. de renormalizacion
!
!
program eje4a
implicit none
  integer, parameter ::d=selected_real_kind(16,200)
  integer :: ii
  real(d) :: probb,proab,proab2,proab4
  real(d) :: norm,norm2,norm4
  real(d) :: ff0,var,pp,suma,suma2,suma4
  real(d), dimension(6) :: dat

  data dat/2._d,1.3_d,3._d,1.5_d,2._d,1.8_d/

  write(*,*)'Ingrese potencia de prob(B): '
  read(*,*)pp
  var=1._d

! calcula las ctes. de renormalizacion
!
  ff0=1._d
  norm=0._d
  norm2=0._d
  norm4=0._d
  do while(ff0 < 4._d)
    probb=ff0**(pp)
    suma=0._d
    do ii=1,6
      suma=suma+(dat(ii)-ff0)**2
      if(ii == 2)suma2=suma
      if(ii == 4)suma4=suma
    enddo
    norm=norm+probb*exp(-suma/2._d/var)
    norm2=norm2+probb*exp(-suma2/2._d/var)
    norm4=norm4+probb*exp(-suma4/2._d/var)
    
    ff0=ff0+0.05_d
  enddo

! ahora hace el calculo
!
  open(7,file='eje4a.out')
  ff0=1._d
  do while(ff0 < 4._d)
    probb=ff0**(pp)
    suma=0._d
    do ii=1,6
      suma=suma+(dat(ii)-ff0)**2
      if(ii == 2)suma2=suma
      if(ii == 4)suma4=suma
    enddo
    proab=probb*exp(-suma/2._d/var)
    proab2=probb*exp(-suma2/2._d/var)
    proab4=probb*exp(-suma4/2._d/var)
    write(7,*)ff0,proab2/norm2,proab4/norm4,proab/norm
    ff0=ff0+0.05_d
  enddo
  close(7)
  
end program eje4a
