Rigorous Computing of Rectangle Scan Probabilities for Markov Increments
Extending recent work of Corrado, we derive an algorithm that computes rigorous upper and lower bounds for rectangle scan probabilities for Markov increments. We experimentally examine the closeness of the bounds computed by the algorithm and we exam…
Authors: Jannis Dimitriadis
RIGOROU S COMPUTING OF R E CT ANGLE SCAN PR OBABILITIE S F OR MAR KO V INCREME N TS By Jannis Dimitriadis Universit¨ at T ri e r April 18, 2019 Rectan gle.Sc an.Probabilities.20110914.tex Extending rec en t work of Cor rado, we deriv e an a lgorithm that computes r igorous uppe r and lower bounds fo r rectangle scan probabilities for Markov increments. W e expe riment ally examine the closeness of the bo unds co mputed by the algor ithm and we exa mine the range of trac table input v aria bles. 1. Introduct ion. Let n balls ra ndomly fall into d b oxe s, each ball w ith probabilit y p i in to b ox i ∈ { 1 , . . . , d } , indep enden tly f rom all the other balls. What is the probability that there exist ℓ adjacen t b o xes in whic h together lie more than k balls? F or ma lly , if w e turn to compute the probability of the complemen t: Let N ∼ M n,p b e a multinomially distributed random v aria ble. T ask: Compute P ( N 1 + . . . + N ℓ ≤ k , . . . , N d − ℓ +1 + . . . + N d ≤ k ) In this pap er, we derive an algorithm that allows fa st computat ion of this probabilit y . Suc h probabilities are needed as p-v alues for tests that c hec k data on clus- ters. F or example: Let n = 50 0 patien t s arrive a t a clinic in d = 365 da ys. W e compute t he probability that there exist three successiv e da ys in whic h together mo r e than 15 patients ar r iv e. F rom the line for k = 15 in T able 3 on page 11 b elo w, w e get the approximate v alue 1 − 0 . 9979961 = 0 . 002003 9 with an a bsolute error less than 10 − 7 . As this probabilit y is so small w e w o uld, if the described ev ent o ccurs, reject the hy p othesis that the patien ts arriv ed indep enden tly and hence susp ect that there m ust b e a reason fo r t his cluster. The supp ort D = { x ∈ N d 0 : x 1 + . . . + x d = n } of the m ultinomial distribution M n,p is finite. Hence w e could compute the desired probabilit y as follow s: F or eac h x ∈ D with x 1 + . . . + x ℓ ≤ k , . . . , x d − ℓ +1 + . . . + x d ≤ k compute the probabilit y P ( N = x ) = n ! / ( x 1 ! . . . x d !) p x 1 1 . . . p x d d and sum up AMS 2000 subje ct classific ations: P rimary 62 P10. Keywor ds and phr ases: Scan Statistics, Recta ngle Probabilities , Marko v Chain, Interv al Arithmetic. 1 2 J. DIMITRIADIS, APRIL 18, 2 019 these v alues. But because the supp ort D is large, this pr o cedure tak es muc h time. F or example : If n = 15 , d = 12 it to ok 4 seconds to compute the probabilit y π := P ( N 1 + N 2 + N 3 ≤ 5 , . . . , N d − 2 + N d − 1 + N d ≤ 5) on a 3 GHz desktop pc, for n = 15 , d = 25 it already to ok 8-9 hours. T o deriv e a f aster metho d in Sections 2 and 3 of this pap er, w e use a fact already utilized b y Corrado [3], na mely that the m ultinomially distributed random v ar ia ble N is a Mark o v incremen t, see Section 4. In this pap er, a Mark ov incr ement is a vec tor ( Y 1 , . . . , Y d ) of discretely distributed ran- dom v ariables with v alues in a group ( X , · ) with the prop ert y that ( Y 1 , Y 1 · Y 2 , . . . , Y 1 · · · Y d ) is a Mark o v c ha in. Our method a ctually w orks for Mark o v incremen ts in this generalit y . F or example, the computation of the probability π for n = 15 , d = 2 5 with the new metho d tak es less than one second. In Sections 5 and 6 w e turn to computer-implemen tations of our alg orithm within the IEEE-754-standard [2] for floating p oint computer arithmetic. The floating p oint n um b er systems according t o the IEEE-754-standard that usual computers w o rk with ha v e the follo wing prop erties: The exact result of an op eration on t w o floating p oin t num b ers, e.g. addition, need not b e a float- ing p oin t n um b er again. In that case , the computer returns a floating p oin t n um b er that is as close as p ossible to the exact result. The difference b et w een the returned v alue and t he exact result is called rounding error . Because of rounding errors, computed v alues, e.g. probabilities, are usually just ap- pro ximations for the exact v alues and the go o dness of the appro ximation is not kno wn. One can switc h the r ounding mo de of the mac hine in suc h a w ay , that in ev ery op eration it returns the minimal floating p oint num b er which is greater or equal than the exact result. This “rounding up” mo de can b e used to compute upp er b ounds fo r the exact v alue, if only p ositiv e num b ers o ccur and only additions and m ultiplications are perfo rmed. In the same w ay , per “rounding do wn” mo de, lo w er b ounds can b e computed. Th us one gets an in- terv al whose b o unds a r e floating p o in t num b ers and in whic h the exact v alue is know n to lie. The accuracy of the a pproximations can easily b e estimated, b ecause the tw o b ounds of the interv al are kno wn. In Section 6, w e presen t an implemen tation of our algorithm within R . F or definiteness, w e assume that all computations are do ne in double-precision according to the IEEE- 754-standard. W e analyze the accuracy of t he R -implemen ta tion and compare it to the b est p ossible accuracy in IEEE-Double-Precision-computations of probabilities, whic h w e examine in Section 5. T o sum up: This w ork extends Corrado‘s by clarifying the underlying Mark ov incremen t structure, b y allowing the computation o f scan proba- bilities and b y providing rigoro us nume rical b o unds. RECT ANGLE SCAN PROBABILITIES FOR MARK OV INCREME NTS 3 2. An algorithm that computes rec t angle probabilities for Marko v incremen t s. W e deriv e an algorithm that computes rectangle probabilities for Marko v incremen ts. It is based on the follo wing recursion form ula: Theorem 2.1 . L et Y = ( Y k ) d k =1 b e Markov in c r ement of a Markov chain ( X k ) d k =1 which takes values in a gr oup ( X , · ) . L et A 1 , . . . , A d ⊂ X b e c ountable sets. Then the pr ob a b i l i ties p ( k , x ) := P ( X k = x, Y 1 ∈ A 1 , . . . , Y k ∈ A k ) for k ∈ { 1 , . . . , d } and x ∈ X fulfil l the r e cursion (1) p ( k , x ) = X y ∈ A k P ( X k = x | X k − 1 = xy − 1 ) p ( k − 1 , xy − 1 ) for k ≥ 2 . Her e and thr oughout, we use the c onv e ntion P ( A | B ) = P ( A ∩ B ) / P ( B ) := 0 if P ( B ) = 0 . Proo f. The functions f k : X 2 → X define d by f k ( x 1 , x 2 ) = x − 1 1 x 2 ha v e the prop ert y that Y k = f k ( X k − 1 , X k ) and f k ( · , x ) is bijectiv e for ev ery x ∈ X . Using this (whic h is actually all w e need, so the metho d w orks not only for Mark ov incremen ts but actually fo r any f unctions of tw o successiv e states o f a Mark o v c hain having the ab ov e bijectivit y prop ert y) a nd writing g k ( x, · ) := f k ( · , x ) − 1 , w e get: P ( X k = x, Y 1 ∈ A 1 , . . . , Y k ∈ A k ) = X y ∈ A k P ( X k = x, Y k = y , Y 1 ∈ A 1 , . . . , Y k − 1 ∈ A k − 1 ) = X y ∈ A k P ( X k = x, X k − 1 = g k ( x, y ) , Y 1 ∈ A 1 , . . . , Y k − 1 ∈ A k − 1 ) = X y ∈ A k P ( X k = x | X k − 1 = g k ( x, y )) × P ( X k − 1 = g k ( x, y ) , Y 1 ∈ A 1 , . . . , Y k − 1 ∈ A k − 1 ) In t he last step the Mark ov prop ert y w as use d. F rom the recursion form ula we c an deriv e the fo llo wing a lgorithm that computes the probability P ( Y 1 ∈ A 1 , . . . , Y d ∈ A d ). Let A 1 , . . . , A d b e finite, so that w e get a finite algorithm. 4 J. DIMITRIADIS, APRIL 18, 2 019 Algorithm A: 1. F or ev ery x ∈ A 1 compute the v alue p (1 , x ) = P ( X 1 = x ) 2. F or ev ery k ∈ { 2 , . . . , d } : F or ev ery x ∈ A 1 · . . . · A k compute the v alue p ( k , x ) with form ula (1) 3. Compute P ( Y 1 ∈ A 1 , . . . , Y d ∈ A d ) = X x ∈ A 1 · ... · A d P ( X d = x, Y 1 ∈ A 1 , . . . , Y d ∈ A d ) Here, let A 1 · · · A n := { a 1 · · · a n : a 1 ∈ A 1 , . . . , a n ∈ A n } , if X is a group and A 1 , . . . , A n ⊂ X . 3. Computing rectangle scan probabilities for Mark ov incremen ts. In t his section w e describe ho w to compute a rectangle s can p robabilit y q := P ( Y 1 · . . . · Y ℓ ∈ A 1 , . . . , Y d − ℓ +1 · . . . · Y d ∈ A d − ℓ +1 ) for a Mark ov incremen t Y . W e us e the following obv ious and w ell- known lemm a: Lemma 3.1 . L et X b e a c ountable set an d ( X k ) d k =1 an X -value d Markov chain. L et W k := ( X k , . . . , X k + ℓ − 1 ) . Then ( W k ) d − ℓ +1 k =1 is an X ℓ -value d Markov chain with tr ansition pr ob abi l i ties P ( W k +1 = w | W k = v ) = P ( X k + ℓ = w ℓ | X k + ℓ − 1 = v ℓ ) for v , w ∈ X ℓ with P ( W k = v ) > 0 and v 2 = w 1 , . . . , v ℓ = w ℓ − 1 . The desired rectangle scan probability for the Mark o v incremen t Y can b e written as a rectangle probabilit y for the incremen t V o f W : If w e set B k := { ( y 1 , . . . , y ℓ ) ∈ X ℓ | y 1 · . . . · y ℓ ∈ A k } we hav e q = P ( V 1 ∈ B 1 , . . . , V d − ℓ +1 ∈ B d − ℓ +1 ) b ecause V k = ( X k − X k − 1 , . . . , X k + ℓ − 1 − X k + ℓ ) for k ∈ { 2 , . . . , d − ℓ + 1 } . The s ets B 1 , . . . , B d − ℓ +1 are p ossibly infinite so t he Algorithm A f rom the last section w ould not work. But if there exist finite sets M 1 , . . . , M d − ℓ +1 ⊂ X ℓ with P ( V 1 ∈ B 1 , . . . , V d − ℓ +1 ∈ B d − ℓ +1 ) = P ( V 1 ∈ M 1 , . . . , V d − ℓ +1 ∈ M d − ℓ +1 ) w e can apply the Algorithm A and th us a re able to c ompute the desired probabilit y . RECT ANGLE SCAN PROBABILITIES FOR MARK OV INCREME NTS 5 Example: If X = ( Z , +) and Y is a Mark o v incremen t with Y 1 , . . . , Y d ≥ 0, then f o r finite sets A 1 , . . . , A d − ℓ +1 ⊂ Z the probability P ( Y 1 + . . . + Y ℓ ∈ A 1 , . . . , Y d − ℓ +1 + . . . + Y d ∈ A d ) equals P (( Y 1 , . . . , Y ℓ ) ∈ M 1 , . . . , ( Y d − ℓ +1 , . . . , Y d ) ∈ M d − ℓ +1 ) with M k := { ( y 1 , . . . , y ℓ ) ∈ N ℓ 0 | y 1 + . . . + y ℓ ∈ A k } , whic h are finite. 4. Examples for Mark ov incremen ts: Multinomially and m ulti- v ariate h ypergeometr ically distributed random v ectors. By b n,p ( k ) = n k p k (1 − p ) n − k w e denote the binomial dens it y with pa r ameters n ∈ N and p ∈ [0 , 1]. By h n,r,b ( k ) = r k b n − k / r + b n w e denote the hy p ergeometrical den- sit y with parameters r , b ∈ N 0 and n ∈ { 1 , . . . , r + b } . Multinomially distributed rando m v ectors as well a s multiv ariate h yp er- geometrically distributed random v ectors are Mark ov increm en t s, hence the results from the last t w o sections a r e applicable in these cases. More precisely , w e ha v e the follo wing tw o propo sitions, as easy calculation with densit y for- m ulas a nd cancelling yield. Example 4.1 . L et ( N 1 , . . . , N d ) ∼ M n,p b e a multinomial ly distribute d r andom varia ble and S k := P k i =1 N i . Then ( S 1 , . . . , S d ) is a Markov chain with P ( S k +1 = x | S k = y ) = b n − y ,p k +1 / P d i = k +1 p i ( x − y ) Example 4.2 . L et ( N 1 , . . . , N d ) ∼ H n, ( m 1 ,...,m d ) b e a multivariate hyp er- ge ometric al ly distribute d r an dom variable, i.e . P ( N 1 = k 1 , . . . , N d = k d ) = m 1 k 1 . . . m d k d / m 1 + ... + m d n for k 1 ∈ { 0 , . . . , m 1 } , . . . , k d ∈ { 0 , . . . , m d } with k 1 + . . . + k d = n , and S k := P k i =1 N i . Then ( S 1 , . . . , S d ) is a Markov chain with P ( S k +1 = x | S k = y ) = h n − y ,m k , P d i = k +1 m i ( x − y ) 5. Definitions and n otations for accuracy analyses of algorithms. In this section w e define terms we need to precisely describe the b eha viour and the accuracy o f n umerical algorithms. F or M ⊂ ]0 , ∞ [ let − M := {− x : x ∈ M } and ± M := M ∪ ( − M ). The IE EE-Double-Precision-Num b er-System is the set IEEE-Double := ± F ∪ ± G ∪ { 0 , −∞ , ∞} 6 J. DIMITRIADIS, APRIL 18, 2 019 with F := { m · 2 e : m ∈ { 2 52 , . . . , 2 53 − 1 } , e ∈ {− 1074 , . . . , 971 }} and G := { k · 2 − 1074 : k ∈ { 1 , . . . , 2 52 − 1 }} , compare [2]. The v alues k / 2 52 in the defi- nition of G and the v alues ( m − 2 52 ) / 2 52 in the definition of F a re called man- tissas of the considered IEEE-Double-Num b ers. W e consider the calculation of probabilities on computation systems that use IEEE-Double-Precision- Num b ers. Hence, eve ry computable pro ba bilit y lies in t he set IEEE-D ouble ∩ [0 , 1] = G ∪ { m · 2 e : m ∈ { 2 52 , . . . , 2 53 − 1 } , e ∈ {− 1074 , . . . , − 53 }} ∪ { 0 , 1 } , the minimal computable probability whic h is greater than zero is min { x ∈ IEEE-Double : x > 0 } = 2 − 1074 ≈ 5 · 10 − 324 and the maximal computable probabilit y wh ic h is les s than one is max { x ∈ IEEE-Double : x < 1 } = 1 − 2 − 53 ≈ 1 − 10 − 16 . W e fix an ob ject not b elonging to the set IEEE-Double, call it NaN for ”Not a Num b er”, and de fine the four operatio ns + , + , · , · : IEEE- Double → IEEE-Do uble ∪ { NaN } F or x, y ∈ IEEE-Double and ◦ ∈ { + , ·} : x ◦ y := min { z ∈ IEEE-Double : z ≥ x ◦ y } x ◦ y := max { z ∈ IEEE-Double : z ≤ x ◦ y } except for the follo wing cases: If x = 0 and y ∈ {−∞ , ∞} or y = 0 and x ∈ { −∞ , ∞} then x · y := x · y := NaN. If x = −∞ and y = ∞ or y = −∞ and x = ∞ then x + y := x + y := NaN. Note t ha t the asso ciativ e law do es not hold for these fo ur op erations. F o r example let a = − 1 , b = 1 , c = 2 − 53 , then we hav e a + ( b + c ) = 0 6 = 2 − 53 = ( a + b )+ c . F or t he calculation of error b o unds for the Algor ithm A deriv ed in Sec- tion 2, w e use the following simple fact: Lemma 5.1 . L et ◦ ∈ { + , ·} , x, y ∈ ]0 , ∞ [ and b 1 , b 2 , c 1 , c 2 ∈ IEEE - Do uble with b 1 ≤ x ≤ c 1 and b 2 ≤ y ≤ c 2 . Then b 1 ◦ b 2 ≤ x ◦ y ≤ c 1 ◦ c 2 F or a quan titativ e ana lysis of the a ccuracy of computed probabilities we need to consider absolute and relativ e errors. F or p, ˜ p ∈ [0 , 1] w e define the absolute error e abs ( p, ˜ p ) := | p − ˜ p | and the relativ e error e rel ( p, ˜ p ) := max e abs ( p, ˜ p ) p , e abs (1 − p, 1 − ˜ p ) 1 − p = | p − ˜ p | min( p, 1 − p ) RECT ANGLE SCAN PROBABILITIES FOR MARK OV INCREME NTS 7 in the approxim ation of p b y ˜ p , with 0 0 := 0 a nd x 0 := ∞ for x > 0. F or a, b ∈ [0 , 1] w ith a ≤ b and ˜ p ∈ [ a, b ] w e f urther define the a bsolute error e abs ([ a, b ] , ˜ p ) := max p ∈ [ a,b ] e abs ( p, ˜ p ) = max { b − ˜ p, ˜ p − a } and the relativ e error e rel ([ a, b ] , ˜ p ) := max p ∈ [ a,b ] e rel ( p, ˜ p ) in the approximation of a probabilit y whic h is kno wn to lie in [ a, b ] b y ˜ p . W e get simple formulas for e rel ([ a, b ] , ˜ p ) in the fo llowing tw o cases. If a, b ∈ [0 , 1 / 2] or a, b ∈ [1 / 2 , 1] we hav e e rel ([ a, b ] , ˜ p ) = max { e rel ( a, ˜ p ) , e rel ( b, ˜ p ) } Hence, if a, b ∈ ]0 , 1 / 2] w e ha ve e rel ([ a, b ] , ˜ p ) = max { ˜ p − a a , b − ˜ p b } and if a, b ∈ [1 / 2 , 1[ w e ha v e e rel ([ a, b ] , ˜ p ) = max { ˜ p − a 1 − a , b − ˜ p 1 − b } F or accuracy measuremen ts in interv al calculations we use the fo llo wing mini- max erro rs: Definition 5.1 . F o r a, b ∈ [0 , 1 ] with a ≤ b we define the absolute error e abs ([ a, b ]) := min ˜ p ∈ [ a,b ] e abs ([ a, b ] , ˜ p ) = e abs ([ a, b ] , a + b 2 ) = b − a 2 and the relativ e error e rel ([ a, b ]) := min ˜ p ∈ [ a,b ] e rel ([ a, b ] , ˜ p ) in the appr oximation of a pr ob ability by the interval [ a, b ] . Easy calculations yield the followin g form ulas: 8 J. DIMITRIADIS, APRIL 18, 2 019 Theorem 5.1 . If a, b ∈ [0 , 1 / 2] we have ∀ ˜ p ∈ [ a, b ] : e rel ([ a, b ] , ˜ p ) ≤ e rel ([ a, b ] , 2 ab a + b ) = b − a b + a Henc e e rel ([ a, b ]) = b − a b + a If a, b ∈ [1 / 2 , 1] we have ∀ ˜ p ∈ [ a, b ] : e rel ([ a, b ] , ˜ p ) ≤ e rel ([ a, b ] , a + b − 2 ab 2 − a − b ) = b − a 2 − a − b Henc e e rel ([ a, b ]) = b − a 2 − a − b Note that the absolute error e abs ([ a, b ]) and the relativ e error e rel ([ a, b ]) need no t b e reache d sim ultaneously by one of the appro ximators. It need not b e reache d at all, as the follo wing example illustrates. Example 5.1 . In T a ble 1 we liste d the err ors e abs ([ a, b ] , ˜ p ) and e rel ([ a, b ] , ˜ p ) for [ a, b ] = [0 . 02 , 0 . 03] and d iffer ent ap pr oximators ˜ p . We se e that e abs ([ a, b ]) = 0 . 005 an d e rel ([ a, b ]) = 1 / 5 . If we take the upp er b ound ˜ p = b as appr oximator for the unknown pr ob ability p , neither e abs ([ a, b ] , ˜ p ) = e abs ([ a, b ]) is r e ache d , nor e rel ([ a, b ] , ˜ p ) = e rel ([ a, b ]) . ˜ p e abs ([ a, b ] , ˜ p ) e rel ([ a, b ] , ˜ p ) 2 ab/ ( a + b ) = 0 . 024 0 . 006 1 / 5 ( a + b ) / 2 = 0 . 0 25 0 . 005 1 / 4 a 0 . 01 1 / 3 b 0 . 01 1 / 2 T able 1 If, fo r ex a mple, the unknown pr ob ability is p = (3 / 10) 3 = 0 . 027 , then the err ors ar e as liste d in T able 2. W e study the maximal accuracy reac hable in double-precision probability calculations: RECT ANGLE SCAN PROBABILITIES FOR MARK OV INCREME NTS 9 ˜ p e abs ( p, ˜ p ) e rel ( p, ˜ p ) 0 . 024 0 . 003 3 / 27 0 . 025 0 . 002 2 / 27 a 0 . 007 7 / 27 b 0 . 003 3 / 27 T able 2 Definition 5.2 . The maximal accuracy in a d ouble-pr e cision c alcula- tion of a pr ob ability p ∈ [0 , 1] is e rel ( p ) := e rel ( I ( p )) wher e I ( p ) := [max { x ∈ IEEE - Double : x ≤ p } , min { x ∈ IEEE - Double : x ≥ p } ] is the minimal interval c ontaining p , w hose endp oi n ts ar e I EEE-Double-Pr e cis i o n-Numb ers. Easy calculation yields e rel ( p ) = ∞ 0 < p < 2 − 1074 or 1 − 2 − 53 < p < 1 1 2 m +1 m ∈ { 1 , . . . , 2 52 − 1 } , p ∈ 2 − 1074 · ] m, m + 1 [ or p ∈ 1 − 2 − 53 · ] m, m + 1 [ 1 2 m +1 m ∈ { 2 52 , . . . , 2 53 − 1 } , e ∈ {− 1022 , . . . , − 2 } , p ∈ 2 e − 52 · ] m, m + 1[ 0 p ∈ IEEE − Double ∩ [0 , 1] F rom the last form ula it follows that • in IEEE-Double floatingp o in t arithmetic w e a re able to approx imate probabilities in [2 − 1074 , 1 − 2 − 53 ] ∪ { 0 , 1 } with finite relative erro r e rel ( p ). • for p ∈ [2 − 1022 , 1 / 2 ], the maximal accuracy satisfies e rel ( p ) ≤ 1 / (2 53 + 1) ≈ 1 . 11 · 10 − 16 . 6. R Implemen tation of Mark o v incremen t scan algorit hms in in- terv al arithmetic. R is an op en source softw are for statistical computa- tions. W e extended R by a C-function that, as p er C-Standard [1] and IEEE- 754-Standard [2], allow s the op erations on IEEE-Double-Num bers whic h w e defined in the previous section. W e wrote an R -progra m t ha t impleme n ts the Algorithm A from Section 2 and us es the principle stated in Lemma 5.1 to compute b o unds for rectangle scan probabilities for Marko v incre- men ts. W e impleme n ted the multinomial a nd m ultiv ariate h yp ergeometric transition pro babilities, as described in Section 4. In a last step the resulting R -implemen tatio n of Algorithm A sets t he returned v alue to 1, if the original return v alue is greater than 1. 6.1. Examples. F or N ∼ M n,p with n = 50 0, d = 365, p = (1 /d , . . . , 1 / d ) and k ∈ { 4 , . . . , 32 } we computed an upp er b ound p and a low er b ound p for 10 J. DIMITRIADIS, APRIL 18, 2 019 the probability P (max d − 2 i =1 ( N i + N i +1 + N i +2 ) ≤ k ) b y an R -implemen ta t io n of the Algo rithm A f rom Section 2. In T able 3 w e tabulate the computed b ounds p , p and analyze their accuracy . Num b ers written in t yp ewriter fon t are hexadecimal. The coloumn titled “ appro x” g iv es the kno wn decimal digits of a v alue o f the “probabilit y repres en ta tion n umber system” T , that lies nearest to the exact v alue. The pro ba bilit y represen tation n um b er s ystem T consists o f all n um b ers with 7 decimal digits without leading z eros or nines. W e use the nota tion . 0 x as an abbreviation for a decimal p oin t follow ed by x zeros, ana lo gously . 9 x . T he sym b ol ? app earing in a n um b er means that the follo wing dig it s are not exactly kno wn. The v a lue e abs resp. e rel is the minimal upp er b ound for e abs ([ p , p ]) resp. e rel ([ p, p ]) which has the form c · 10 k where c has 3 significan t dig it s and k ∈ Z . Th us, the line with k = 15 means that the probability P (max d − 2 i =1 ( N i + N i +1 + N i +2 ) ≤ 15) lies in the in terv al [ p , p ] with p = 1.fef956911fe58 · 2 − 1 = (1 + 15 · 1 6 − 1 + 14 · 16 − 2 + . . . + 8 · 16 − 13 ) · 2 − 1 = 0 . 9979 960491327 330984745 4584087245166301727 2949 2 1875 p = 1.fef95690c7eda · 2 − 1 = (1 + 15 · 1 6 − 1 + . . . + 10 · 16 − 13 ) · 2 − 1 = 0 . 9979 960490927 297644958 5712543921545147895 8129 8 828125 with all equalities exact. The minimal upp er bound for e abs ([ p , p ]) whic h has the form c · 10 k where c has 3 significan t digits and k ∈ Z is 2 . 01 · 10 − 11 and the minimal upp er b o und for e abs ([ p, p ]) whic h has this fo r m is 9 . 9 9 · 10 − 9 . A v alue of the n umber sys tem T whic h is nearest to the ex act probabilit y is 0 . 9979961 . As the n um b ers of the system T in the interv a l [0 . 001 , 0 . 9989999] differ b y 10 − 7 , just kno wing the appro ximate v alue w e can infer that the absolute erro r in this appro ximation is less than 10 − 7 . 6.2. R em arks on numeric al c omputations of multinomial pr ob abi l i ties. 6.2.1. R elative err or of c o m plement pr ob a bilities. In the preceding sec- tion w e computed the distribution function of a multinomial scan statistic. F or sev eral applications, e.g. m ultinomial scan te st, w e w a nt to compute the upp er distribution function ins tead. If w e compute its v alues from t he com- plemen ts P (max d i =1 N i + N i +1 + N i +2 ≥ k ) = 1 − P ( max d i =1 N i + N i +1 + N i +2 ≤ k − 1) in exact arithmetic, for example b y using a suitable softw are, there is no increse of erro r. I f we do auto ma t ic computation of the compleme n t in IEEE-Double-Precision-Num b er- System, then the error increases for small RECT ANGLE SCAN PROBABILITIES FOR MARK OV INCREME NTS 11 k p, p e abs e rel appro x 4 0 0 0 0 5 1.1c5df1 e1a1f83 · 2 − 178 1.1c5df1 e171043 · 2 − 178 5 . 82 · 10 − 65 2 . 01 · 10 − 11 . 0 53 28993 6 1.b826f2 2f10057 · 2 − 67 1.b826f2 2ec43c3 · 2 − 67 2 . 34 · 10 − 31 2 . 01 · 10 − 11 . 0 19 11651 7 1.b71c49 2587c97 · 2 − 27 1.b71c49 253c2df · 2 − 27 2 . 57 · 10 − 19 2 . 01 · 10 − 11 . 0 7 12780 8 1.98b835 1d76fbd · 2 − 11 1.98b835 1d309cf · 2 − 11 1 . 57 · 10 − 14 2 . 01 · 10 − 11 . 0 3 77957 9 1.0f0230 ce6f8a1 · 2 − 4 1.0f0230 ce40e15 · 2 − 4 1 . 33 · 10 − 12 2 . 01 · 10 − 11 . 06616 42 10 1.826e2a db7befd · 2 − 2 1.826e2a db39686 · 2 − 2 7 . 57 · 10 − 12 2 . 01 · 10 − 11 . 37737 34 11 1.7131cf 887a229 · 2 − 1 1.7131cf 883a935 · 2 − 1 1 . 45 · 10 − 11 5 . 19 · 10 − 11 . 72108 32 12 1.ce5760 94ddb84 · 2 − 1 1.ce5760 948e1f6 · 2 − 1 1 . 81 · 10 − 11 1 . 87 · 10 − 10 . 90301 04 13 1.f11623 01d80ec · 2 − 1 1.f11623 01827ae · 2 − 1 1 . 95 · 10 − 11 6 . 69 · 10 − 10 . 97087 20 14 1.fbef94 98b0df9 · 2 − 1 1.fbef94 98596d7 · 2 − 1 1 . 99 · 10 − 11 2 . 51 · 10 − 9 . 99206 22 15 1.fef956 911fe58 · 2 − 1 1.fef956 90c7eda · 2 − 1 2 . 01 · 10 − 11 9 . 99 · 10 − 9 . 99799 61 16 1.ffc1fb bfd6e58 · 2 − 1 1.ffc1fb bf7ecb1 · 2 − 1 2 . 01 · 10 − 11 4 . 24 · 10 − 8 . 9 3 52685 17 1.fff23b 0d23a3c · 2 − 1 1.fff23b 0ccb810 · 2 − 1 2 . 01 · 10 − 11 1 . 91 · 10 − 7 . 9 3 89495 18 1.fffd1d 22cb527 · 2 − 1 1.fffd1d 22732da · 2 − 1 2 . 01 · 10 − 11 9 . 11 · 10 − 7 . 9 4 77980 19 1.ffff6d 5024936 · 2 − 1 1.ffff6d 4fcc6e4 · 2 − 1 2 . 01 · 10 − 11 4 . 59 · 10 − 6 . 9 5 56284 20 1.ffffe4 570f39a · 2 − 1 1.ffffe4 56b7146 · 2 − 1 2 . 01 · 10 − 11 2 . 44 · 10 − 5 . 9 6 17567 21 1.fffffb 08bd13c · 2 − 1 1.fffffb 0864ee9 · 2 − 1 2 . 01 · 10 − 11 1 . 36 · 10 − 4 . 9 6 8520 ? 22 1.ffffff 264f47d · 2 − 1 1.ffffff 25f7228 · 2 − 1 2 . 01 · 10 − 11 7 . 91 · 10 − 4 . 9 7 74 ? 23 1.ffffff dc79315 · 2 − 1 1.ffffff dc210c0 · 2 − 1 2 . 01 · 10 − 11 4 . 83 · 10 − 3 . 9 8 6 ? 24 1.ffffff fa913ba · 2 − 1 1.ffffff fa39167 · 2 − 1 2 . 01 · 10 − 11 3 . 08 · 10 − 2 . 9 9 ? 25 1.ffffff ff53a50 · 2 − 1 1.ffffff fefb7fe · 2 − 1 2 . 01 · 10 − 11 2 . 04 · 10 − 1 . 9 9 ? 26 1 1.ffffff ffb44b7 · 2 − 1 1 − p ∞ . 9 10 ? 27 1 1.ffffff ffcf373 · 2 − 1 1 − p ∞ . 9 10 ? 28 1 1.ffffff ffd2fd3 · 2 − 1 1 − p ∞ . 9 10 ? 29 1 1.ffffff ffd37fa · 2 − 1 1 − p ∞ . 9 10 ? 30 1 1.ffffff ffd3908 · 2 − 1 1 − p ∞ . 9 10 ? 31 1 1.ffffff ffd392a · 2 − 1 1 − p ∞ . 9 10 ? 32 1 1.ffffff ffd392a · 2 − 1 1 − p ∞ . 9 10 ? T able 3 Upp er and lower b oun ds p, p for P (max d − 2 i =1 N i + N i +1 + N i +2 ≤ k ) with N ∼ M n,p , n = 5 00 , d = 365 , p = (1 /d, . . . , 1 /d ) and k ∈ { 4 , . . . , 32 } . F or details, se e Subse ction 6.1. 12 J. DIMITRIADIS, APRIL 18, 2 019 probabilities. Then, we are not able to approximate probabilities less then 10 − 16 with a finite relativ e error. Compare T able 4. Here, the relativ e er- ror increas es for small probabilities as w ell as for big probabilities. Small complemen ts of probabilities are lo st. In general, one should try to a v o id de - v eloping a lg orithms that complemen t t he computed probability at the end. An algor it hm that computes the complemen t is not equiv alen t to a direct one. The maximal accuracy with resp ect to complemen t ation o f proba bilities is defined as e rel , c ( p ) := max( e rel ( p ) , e rel (1 − p )) Easy calculation yields e rel ( p ) = ∞ 0 < p < 2 − 53 or 1 − 2 − 53 < p < 1 1 2 m +1 m ∈ { 1 , . . . , 2 52 − 1 } , p ∈ 2 − 53 · ] m, m + 1[ or p ∈ 1 − 2 − 53 · ] m, m + 1[ 0 p ∈ IEEE − Double ∩ [0 , 1] 6.2.2. Co m putation Time and Sp ac e. Besides the a ccuracy o f the a lgo- rithm, there are tw o other problems t hat matter: Time and space needed to compute the probabilit y . The imple men tat io n of the m ultinomial scan algorithm w e made needs to store 2 ∗ n + ℓ ℓ Double-Precision-Num bers. Eac h Double-Precision-Num b er needs 8 Bytes. F or example, for the scan width ℓ = 3, on a computer with 16 GByte memory , w e were able to compute Scan-Pro ba bilities for up t o ap- pro ximately n = 1700 in double-precision. Using Single-Precision-Num b ers, whic h take o nly 4 Bytes, we could compute up t o n = 215 0, but the accuracy is w orse than in double-precision computations, as T able 5 demonstrates. The IE EE-Single-Precision-Num b er-System is the set IEEE-Single := ± F ∪ ± G ∪ { 0 , −∞ , ∞} with F := { m ∗ 2 e : m ∈ { 2 23 , . . . , 2 24 − 1 } , e ∈ {− 149 , . . . , 104 }} a nd G := { k ∗ 2 − 149 : k ∈ { 1 , . . . , 2 23 − 1 }} , compare [2]. In the third coloumn of T a- ble 5 w e listed the first digits of the computed b ounds p , p in decimal format. The t ime that it tak es to compute a rectangle scan proba bility for a m ulti- nomially distributed random vector in single precision do es not differ m uc h from the time it tak es in double-precis ion, examples are listed in T able 6. 6.3. Bin omial Pr ob abilities. W e use the follow ing alg orithm to compute the multinomial transition probabilities, that ar e binomial. RECT ANGLE SCAN PROBABILITIES FOR MARK OV INCREME NTS 13 k p, p e abs e rel appro x 5 1 0 0 1 6 1 1.ffffff fffffff · 2 − 1 1 − p ∞ . 9 15 ? 7 1 1.ffffff fffffff · 2 − 1 1 − p ∞ . 9 15 ? 8 1.ffffff 9238edc · 2 − 1 1.ffffff 9238edb · 2 − 1 5 . 55 · 10 − 17 4 . 34 · 10 − 9 . 9 7 87220 9 1.ff99d1 f2b8b3e · 2 − 1 1.ff99d1 f2b8a24 · 2 − 1 1 . 57 · 10 − 14 2 . 01 · 10 − 11 . 9 3 22042 10 1.de1fb9 e637e3e · 2 − 1 1.de1fb9 e6320eb · 2 − 1 1 . 33 · 10 − 12 2 . 01 · 10 − 11 . 93383 58 11 1.3ec8ea 92634bd · 2 − 1 1.3ec8ea 9242081 · 2 − 1 7 . 57 · 10 − 12 2 . 01 · 10 − 11 . 62262 66 12 1.1d9c60 ef8ad96 · 2 − 2 1.1d9c60 ef0bbae · 2 − 2 1 . 45 · 10 − 11 5 . 19 · 10 − 11 . 27891 68 13 1.8d44fb 5b8f050 · 2 − 4 1.8d44fb 59123e0 · 2 − 4 1 . 81 · 10 − 11 1 . 87 · 10 − 10 . 09698 96 14 1.dd3b9f cfb0a40 · 2 − 6 1.dd3b9f c4fe280 · 2 − 6 1 . 95 · 10 − 11 6 . 69 · 10 − 10 . 02912 80 15 1.041ad9 e9a4a40 · 2 − 7 1.041ad9 d3c81c0 · 2 − 7 1 . 99 · 10 − 11 2 . 51 · 10 − 9 . 00793 77 16 1.06a96f 3812600 · 2 − 9 1.06a96e e01a800 · 2 − 9 2 . 01 · 10 − 11 9 . 99 · 10 − 9 . 00200 40 17 1.f02204 09a7800 · 2 − 12 1.f02201 48d4000 · 2 − 12 2 . 01 · 10 − 11 4 . 24 · 10 − 8 . 0 3 47315 18 1.b89e66 8fe0000 · 2 − 14 1.b89e5b 8b88000 · 2 − 14 2 . 01 · 10 − 11 1 . 91 · 10 − 7 . 0 3 10505 19 1.716ec6 6930000 · 2 − 16 1.716e9a 56c8000 · 2 − 16 2 . 01 · 10 − 11 9 . 11 · 10 − 7 . 0 4 22020 20 1.256067 2380000 · 2 − 18 1.255fb6 d940000 · 2 − 18 2 . 01 · 10 − 11 4 . 59 · 10 − 6 . 0 5 4371 ? 21 1.ba948e ba00000 · 2 − 21 1.ba8f0c 6600000 · 2 − 21 2 . 01 · 10 − 11 2 . 44 · 10 − 5 . 0 6 824 ? 22 1.3de6c4 5c00000 · 2 − 23 1.3dd0bb 1000000 · 2 − 23 2 . 01 · 10 − 11 1 . 36 · 10 − 4 . 0 6 14 ? 23 1.b411bb 0000000 · 2 − 26 1.b36170 6000000 · 2 − 26 2 . 01 · 10 − 11 7 . 91 · 10 − 4 . 0 7 253 ? 24 1.1ef7a0 0000000 · 2 − 28 1.1c3675 8000000 · 2 − 28 2 . 01 · 10 − 11 4 . 83 · 10 − 2 . 0 8 41 ? 25 1.71ba64 0000000 · 2 − 31 1.5bb118 0000000 · 2 − 31 2 . 01 · 10 − 11 3 . 08 · 10 − 2 . 0 9 6 ? 26 1.048020 0000000 · 2 − 33 1.58b600 0000000 · 2 − 34 2 . 01 · 10 − 11 2 . 04 · 10 − 1 . 0 9 ? T able 4 Upp er and lower b oun ds p, p for P (max d − 2 i =1 N i + N i +1 + N i +2 ≥ k ) with N ∼ M n,p , n = 5 00 , d = 365 , p = (1 /d, . . . , 1 /d ) and k ∈ { 5 , . . . , 26 } . 14 J. DIMITRIADIS, APRIL 18, 2 019 k p, p p, p e abs e rel appro x 4 0 0 0 0 0 5 1.974c00 · 2 − 135 0 . 0 40 3652 ... 0 1 . 83 · 10 − 41 1 . 04 · 10 − 2 . 0 40 ? 6 1.bcc5a4 · 2 − 67 1.b39300 · 2 − 67 . 0 19 1177 ... . 0 19 1152 ... 1 . 22 · 10 − 22 1 . 04 · 10 − 2 . 0 19 11? 7 1.bbb862 · 2 − 27 1.b28b40 · 2 − 27 . 0 7 12913 ... . 0 7 12646 ... 1 . 34 · 10 − 10 1 . 04 · 10 − 2 . 0 7 12? 8 1.9d02a2 · 2 − 11 1.947834 · 2 − 11 . 0 3 78775 ... . 0 3 77146 ... 8 . 15 · 10 − 6 1 . 04 · 10 − 2 . 0 3 7? 9 1.11da84 · 2 − 4 1.0c30d0 · 2 − 4 . 0668587 ... . 0654762 ... 6 . 91 · 10 − 4 1 . 04 · 10 − 2 . 06? 10 1.867cac · 2 − 2 1.7e699a · 2 − 2 . 3813349 ... . 3734497 ... 3 . 94 · 10 − 3 1 . 04 · 10 − 2 . 3? 11 1.7511fc · 2 − 1 1.6d5b2a · 2 − 1 . 7286528 ... . 7135861 ... 7 . 53 · 10 − 3 2 . 7 · 10 − 2 . 7? 12 1.d331e6 · 2 − 1 1.c988cc · 2 − 1 . 9124900 ... . 8936218 ... 9 . 43 · 10 − 3 9 . 7 · 10 − 2 . ? 13 1.f64e04 · 2 − 1 1.ebeb16 · 2 − 1 . 9810639 ... . 9607779 ... 1 . 01 · 10 − 2 3 . 49 · 10 − 1 . 9? 14 1 1.f6a7a6 · 2 − 1 1 . 9817478 ... 1 − p ∞ . 9? 15 1 1.f9a956 · 2 − 1 1 . 9876200 ... 1 − p ∞ . 9? 16 1 1.fa6fe6 · 2 − 1 1 . 9891349 ... 1 − p ∞ . 9? 17 1 1.fa9fa0 · 2 − 1 1 . 9894990 ... 1 − p ∞ . 9? 18 1 1.faaa68 · 2 − 1 1 . 9895813 ... 1 − p ∞ . 9? 19 1 1.faacb6 · 2 − 1 1 . 9895989 ... 1 − p ∞ . 9? 20 1 1.faad2c · 2 − 1 1 . 9896024 ... 1 − p ∞ . 9? 21 1 1.faad3c · 2 − 1 1 . 9896029 ... 1 − p ∞ . 9? 22 1 1.faad40 · 2 − 1 1 . 9896030 ... 1 − p ∞ . 9? 23 1 1.faad44 · 2 − 1 1 . 9896031 ... 1 − p ∞ . 9? 24 1 1.faad46 · 2 − 1 1 . 9896032 ... 1 − p ∞ . 9? 25 1 1.faad46 · 2 − 1 1 . 9896032 ... 1 − p ∞ . 9? T able 5 Upp er and lower b oun ds p, p for P (max d − 2 i =1 N i + N i +1 + N i +2 ≤ k ) with N ∼ M n,p , n = 500 , d = 365 , p = (1 /d, . . . , 1 /d ) and k ∈ { 4 , . . . , 25 } , c ompute d in single-pr e cision. RECT ANGLE SCAN PROBABILITIES FOR MARK OV INCREME NTS 15 n d k p (do uble) time (double) p (s ingle) time (sing le) 100 365 6 0 . 9934 578 1 s 0 . 991 4927 1 s 500 365 13 0 . 97087 20 1 min 2 7 s 0 . 960777 9 1 min 2 5 s 1000 3 65 20 0 . 9 60432 4 49 min 39 s 0 . 94055 73 1 0 min 19 s 1500 365 27 0 . 97393 03 3 h 22 min 26 s 0 . 9438 554 2 h 44 min 00 s 1700 3 65 29 0 . 9 61031 5 5 h 47 min 01 s 0 . 92748 42 4 h 58 min 21 s 1750 3 65 31 x no co mputation p ossible 0 . 9516 879 6 h 3 6 min 48 s 2150 365 37 x no computation p oss ible 0 . 95072 57 12 h 47 min 22 s T able 6 Computation time for a lower b oun d p for the s c an pr ob ability P (max d − 2 i =1 N i + N i +1 + N i +2 ≤ k ) for a r andom ve ctor ( N 1 , . . . , N d ) with the mu ltinomial distribution M n,p with p = (1 /d, . . . , 1 /d ) , in s ingle pr e cision and in double pr e cision. Details ar e describ e d in Subsubse ction 6.2.2 double bnp(unsigned int k,unsigned int n, doub le p, double q){ if (2*k>n) return(bnp(n-k ,n,q,p)); double f=1.0; unsigned int j0=0,j1=0,j2 =0; while ( (j0
Original Paper
Loading high-quality paper...
Comments & Academic Discussion
Loading comments...
Leave a Comment