;-----------------------------------------------------------
;idlharj2_poisson_rnd
;HS 050308
;-----------------------------------------------------------

;poisson jakauma:
  lambda=5
  x=findgen(21)
  px=exp(-lambda)*lambda^x/factorial(x)

  nwin
  plot,x,px,xtitle='x',ytitle='p(x)',$
       title='Poisson with lambda='+string(lambda,'(i2)')

;random number using IDL
  r=randomu(seed,100000,poisson=lambda)
  histo_f,r,-0.5,21,1,xx,yy
  oplot,xx,yy,psym=10


;directly
;Poisson process with lambda events/timeunit
;waiting time between successive events is exponentially distributed,
;with a mean waiting time 1/lambda
;add up X+1 exponentially distributed waiting times, until the sum
;exceed timeunit -> X is the number of events during timeunit
; X is Poisson distributed.

  rr=fltarr(100000)
  for i=0l,n_elements(rr)-1 do begin
     z=-1./lambda*alog(randomu(seed,20))
     zc=total(z,/cum)        ;cumulative sum
     ind=where(zc gt 1)
     rr(i)= ind(0)           ;number of event during the time interval
                             ;note first  zc > 1   -> ind(0)=0 (x+1=1, x=0)
                              ;    second zc > 2   -> ind(0)=1 (x+1=2, x=1)
                              ;     etc   zx > x+1 -> ind(0)=x
 endfor
  histo_f,rr,-0.5,21,1,xx,yy
  oplot,xx,yy,psym=10,col=2
  
  
;-------------------------------------------------------------
;illustration
;show 11 trials, each event until TIME> UNIT TIME
;number of event poisson distributed

  nwin
  plot,lindgen(2),/nod,xr=[-.1,3],xs=1,yr=[-1,11],ys=1,xtitle='TIME',$
       ytitle='TRIAL #',$
       title='Poisson process with lambda='+string(lambda,'(i2)')

  for i=0l,10 do begin
     z=-1./lambda*alog(randomu(seed,100))
     zc=total(z,/cum)           ;cumulative sum
     
     ind=where(zc lt 1,count)
     if(count gt 0) then ind=lindgen(count)
     if(count gt 0) then oplot,zc(ind),zc(ind)*0+i,psym=6,col=2,syms=2 
     
     
     oplot,zc(count)*[1,1],[i,i],psym=1,col=3
     
     add=''
     if(i eq 10) then add=' EVENTS/UNIT TIME'
     
     xyouts,zc(count)+.1,i,string(count,'(i2)')+add
     oplot,[0,zc(count)],[i,i],lines=1,col=3
     if(count gt 0) then oplot,[0,zc(count-1)],[i,i],lines=1

   
  endfor
  
;insert of waiting time distribution
  x=findgen(500)/100.*1./lambda
  fx=1./lambda*exp(-x*lambda)
  plot,x,fx,pos=[0.65,0.5,0.9,0.8],/noerase,xtitle='waiting time',ytitle='p'
  oplot,[1.,1.]/lambda,[0,1],lines=2
  
end