gdl-astrolib 2018.02.16+dfsg-1 / usr / share / gnudatalanguage / astrolib / poidev.pro

This file is indexed.

/usr/share/gnudatalanguage/astrolib/poidev.pro is in gdl-astrolib 2018.02.16+dfsg-1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below. 

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
function poidev, xm, SEED = seed
;+
; NAME:
;     POIDEV
; PURPOSE:
;     Generate a Poisson random deviate
; EXPLANATION:
;     Return an integer random deviate drawn from a Poisson distribution with
;     a specified mean.    Adapted from procedure of the same name in 
;     "Numerical Recipes" by Press et al. (1992), Section 7.3
;
;     NOTE: This routine became partially obsolete in V5.0 with the 
;     introduction of the POISSON keyword to the intrinsic functions 
;     RANDOMU and RANDOMN.     However, POIDEV is still useful for adding 
;     Poisson noise to an existing image array, for which the coding is much 
;     simpler than it would be using RANDOMU (see example 1) 
; CALLING SEQUENCE:
;     result = POIDEV( xm, [ SEED = ] )
;
; INPUTS:
;     xm - numeric scalar, vector or array, specifying the mean(s) of the 
;          Poisson distribution
;
; OUTPUT:
;     result - Long integer scalar or vector, same size as xm
;
; OPTIONAL KEYWORD INPUT-OUTPUT:
;     SEED -  Scalar to be used as the seed for the random distribution.  
;             For best results, SEED should be a large (>100) integer.
;             If SEED is undefined, then its value is taken from the system 
;             clock (see RANDOMU).    The value of SEED is always updated 
;             upon output.   This keyword can be used to have POIDEV give 
;             identical results on consecutive runs.     
;
; EXAMPLE:
;     (1) Add Poisson noise to an integral image array, im
;              IDL> imnoise = POIDEV( im)
;
;     (2) Verify the expected mean  and sigma for an input value of 81
;              IDL> p = POIDEV( intarr(10000) + 81)   ;Test for 10,000 points
;              IDL> print,mean(p),sigma(p)
;     Mean and sigma of the 10000 points should be close to 81 and 9
;
; METHOD: 
;     For small values (< 20) independent exponential deviates are generated 
;     until their sum exceeds the specified mean, the number of events 
;     required is returned as the Poisson deviate.   For large (> 20) values,
;     uniform random variates are compared with a Lorentzian distribution 
;     function.
;
; NOTES:
;     Negative values in the input array will be returned as zeros.  
;
;       
; REVISION HISTORY:
;      Version 1               Wayne Landsman        July  1992
;      Added SEED keyword                            September 1992
;      Call intrinsic LNGAMMA function               November 1994
;      Converted to IDL V5.0   W. Landsman   September 1997
;      Use COMPLEMENT keyword to WHERE()        W. Landsman August 2008
;-
  On_error,2
  compile_opt idl2

 Npts = N_elements( xm)
 
 case NPTS of 
 0: message,'ERROR - Poisson mean vector (first parameter) is undefined'
 1: output = lonarr(1) 
 else: output = make_array( SIZE = size(xm), /NOZERO ) 
 endcase 
 
   index = where( xm LE 20, Nindex, complement=big, Ncomplement=Nbig)

   if Nindex GT 0 then begin

   g = exp( -xm[ index] )           ;To compare with exponential distribution
   em1 = replicate( -1, Nindex )    ;Counts number of events
   t = replicate( 1., Nindex )          ;Counts (log) of total time

  Ngood = Nindex
  good = lindgen( Nindex)                 ;GOOD indexes the original array
  good1 = good                         ;GOOD1 indexes the GOOD vector

 REJECT:  em1[good] = em1[good] + 1      ;Increment event counter
   t = t[good1]*randomu( seed, Ngood )   ;Add exponential deviate, equivalent
                                         ;to multiplying random deviate
   good1 = where( t GT g[good], Ngood1)  ;Has sum of exponential deviates 
                                         ;exceeded specified mean?
   if ( Ngood1 GE 1 ) then begin
           good = good[ good1]
           Ngood = Ngood1
           goto, REJECT
   endif
   output[index] = em1
 endif
     if Nindex EQ Npts then return, output
; ***************************************

    xbig = xm[big]

    sq = sqrt( 2.*xbig )           ;Sq, Alxm, and g are precomputed
    alxm = alog( xbig )
    g = xbig * alxm - lngamma( xbig + 1.)

    Ngood = Nbig  & Ngood1 = Nbig
    good = lindgen( Ngood)
    good1 = good
    y = fltarr(Ngood, /NOZERO ) & em = y


REJECT1:   y[good] = tan( !PI * randomu( seed, Ngood ) )  
   em[good] = sq[good]*y[good] + xbig[good]
   good2 = where( em[good] LT 0. , Ngood )
   if (Ngood GT 0) then begin
            good = good[good2]
            goto, REJECT1
   endif

   fixem = long( em[good1] )
   test = check_math( 0, 1)         ;Don't want overflow messages
   t = 0.9*(1. + y[good1]^2)*exp( fixem*alxm[good1] - $ 
               lngamma( fixem + 1.) - g[good1] )
   good2 = where( randomu (seed, Ngood1) GT T , Ngood)
   if ( Ngood GT 0 ) then begin
            good1 = good1[good2]
            good = good1
            goto, REJECT1
   endif
   output[ big ] = long(em)

 return, output

 end

APT Browse - Built by Thomas Orozco.