;--------------------------------------------------
;STATISTICAL METHODS IN ASTRONOMY
; EXERCISES I
; 21.02.2008 Heikki Salo
; exercise1_bino1.pro
;--------------------------------------------------

 program='exercise1_bino1'

;--------------------------------------------------
;calculate binomial probability prob(V;N,p)
;prob(V; N,p) = C_V p^V (1-p)^(N-v)
;
;           N!          N*(N-1)* ... *(N-V+1) 
; C_V =  ---------   =  --------------------
;        V! (N-V)!      1*2* ....       * V

; the latter holds only when V>1   
; C_0 = 1
;--------------------------------------------------

;--------------------------------------------------
;1. with a direct evaluation of binomial coefficients
;           N!       
; C_V =  ---------   
;        V! (N-V)!   
;--------------------------------------------------
;
;--------------------------------------------------
;2. with binomial coefficients calculated recursively
;
; C_k = 1  for k=0
; C_(k+1)=C_k * (N-k)/(k+1)   k=1,2,..N-1
;
;--------------------------------------------------

;--------------------------------------------------
;3. Using the gamma-function
; log(C_v) = LNGAMMA(n+1)-LNGAMMA(v+1)-LNGAMMA(n-v+1)

;--------------------------------------------------

;--------------------------------------------------
;4. Approximating with a Gaussian functian
; Large N --> approaches Gaussian with 
; mu= p*N
; sigma=sqrt(p*(1-p)*N)
; use gauss_pdf which returns cumulative density function of
; normalized gaussian (zero mean, unit standard deviation)

; subsitute  x1= (V-0.5- mu)/sigma
;            x2= (V+0.5- mu)/sigma
; prob(V)  = phi(x2)-phi(x1)   approximatively
;--------------------------------------------------

;-------------------------------------------------
;0. check by bino_pdf  (utilizing idl-library routine)
;-------------------------------------------------
;output to screen (or file)
 ps=0
 psdirect,program,ps,/color
;-------------------------------------------------

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

 p=0.5d0
 N=10

;---------------------------------------------------------
;0. IDL-library routine
  
  vtab0=dindgen(n+1)
  bino_pdf,n,p,vtab0,ptab0

;-----------------------------------------------------------
;1. directly
;calculate and store factorial
;-----------------------------------------------------------
                
  vtab=dindgen(n+1)             ;(largest possible  n about 170)   
  ptab1=vtab*0.

  fac=dindgen(n+1)
  fac(0)=1.
  for i=1,n do begin
     fac(i)=fac(i-1)*i
  endfor 
  for v=0,n do begin
     c_v=fac(n)/fac(v)/fac(n-v)
     ptab1(v)=c_v*p^v*(1.-p)^(n-v)
  endfor

  nwin
  add='   N='+string(n,'(i4)')+' p='+string(p,'(f4.2)')
  plot,vtab,ptab1,title='C_v directly'+ add,psym=6,syms=.5
  oplot,vtab0,ptab0,col=2,lines=2
  
;-----------------------------------------------------------
;2. calculate c_v recursively
;-----------------------------------------------------------

  ptab2=vtab*0.
  c_v=1.d0
  for v=0,n do begin
     ptab2(v)=c_v*p^v*(1.-p)^(n-v)
;binomial factor for the NEXT v
     c_v=c_v*(n-v)/(v+1)
  endfor
  
  nwin
  add='   N='+string(n,'(i4)')+' p='+string(p,'(f4.2)')
  plot,vtab,ptab2,title='C_v recursively'+add,psym=6,syms=.5
  oplot,vtab0,ptab0,col=2,lines=2

;-----------------------------------------------------------
;3. Use gamma-function
;-----------------------------------------------------------

  ptab3=vtab*0.
  for v=0,n do begin     
;remember to make in double precision!
     logfac = LNGAMMA(n+1.d0)-LNGAMMA(v+1.d0)-LNGAMMA(n-v+1.d0)
     ptab3(v)=exp(logfac)*p^v*(1.-p)^(n-v)
  endfor
    
  nwin
  add='   N='+string(n,'(i4)')+' p='+string(p,'(f4.2)')
  plot,vtab,ptab3,title='C_v from gamma-function'+add,psym=6,syms=.5
  oplot,vtab0,ptab0,col=2,lines=2
  
;-----------------------------------------------------------
;4. Use Gaussian approximation
;-----------------------------------------------------------

  ptab4=vtab*0.
  mu=p*n
  sigma=sqrt(p*(1.d0-p)*n)
  
  for v=0,n do begin
     x1= (V-0.5- mu)/sigma
     x2= (V+0.5- mu)/sigma
     ptab4(v)=gauss_pdf(x2)-gauss_pdf(x1)
  endfor
  
  nwin
  add='   N='+string(n,'(i4)')+' p='+string(p,'(f4.2)')
  plot,vtab,ptab4,title='Gaussian approximation'+add,psym=6,syms=.5
  oplot,vtab0,ptab0,col=2,lines=2

 !p.multi=0
  psdirect,program,ps,/stop 
end