5 people like it.

Distribution of Random hyperharmonic series

The random hyperharmonic series is the infinite series S = Sum(1,inf,d(i)/i^pow), where integer pow is greater than 1, and d(i) are independent, identically distributed random variables with property P(d(i)=0) = P(d(i)=1) = 0.5. Cumulative function F(x) = P(S < x) for even powers can be build by combination of analytical and numerical computations.

  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: 
135: 
136: 
137: 
138: 
139: 
140: 
141: 
142: 
143: 
144: 
145: 
146: 
147: 
148: 
149: 
150: 
151: 
152: 
153: 
154: 
155: 
156: 
157: 
158: 
159: 
160: 
161: 
162: 
163: 
164: 
165: 
166: 
167: 
168: 
169: 
170: 
171: 
172: 
173: 
174: 
175: 
(*The random hyperharmonic series is the infinite series S = Sum(1,inf,d(i)/i^pow),
where integer pow > 1, and d(i) are independent, identically distributed 
random variables with property P(d(i)=0) = P(d(i)=1) = 0.5. The fact that 
for even pows the hyperharmonic series converge to some value at [0, ΞΆ(pow)] with prob 1.
To build cumulative function F(x) = P(S < x) the tail of a series replaced by an
appropriate normal random variable. Complexity ~ O(2^k), although faster calculations 
are possible, the program shows connections between number theory and probability theory.
Expressions employed here are quite simple, therefore many improvements can be made easily*)

// some primitive actions
let inline square x = x * x
let inline double x = x + x
let inline half x = x*0.5
let inline inverse x = 1.0/x
let sign x = if x=0.0 then x else x/(abs x)

//Binomial coeffs
type binom =
     static member fact n = 
            let rec loop acc = function
                | n when n = 0 -> acc
                | n -> loop (n*acc) (n-1)
            loop 1 n
   
     static member comb n m = 
            binom.fact n / (binom.fact m * binom.fact (n-m))

//Bernoulli numbers
let rec bernoulli = function
    | n when n = 0 -> 1.0
    | n when n = 1 -> -0.5
    | n when n % 2 = 1 -> 0.0
    | n -> -([0..(n-1)] 
                  |> List.map (fun k -> 
                          (binom.comb (n+1) k |> float, bernoulli k))
                  |> List.sumBy (fun x -> fst x * snd x)   
            ) / (float (n+1))

//optional upper limit for a sum                
type UpperLimit =
   | Infinity
   | Integer of int

//Hyperharmonic series of even power   
type HypHarmonic(pow) =
   let m = pow/2
   let repl = 
       let rec loop f n =
           if n = 1 then f
           else loop (f>>f) (n/2)
       loop square (pow/2) 
 
   member h.pow = pow
   //terms
   member h.a i = i |> repl |> inverse
 
   static member (*) (a:HypHarmonic,b:HypHarmonic) = HypHarmonic(a.pow * b.pow)
   static member pi2 = double System.Math.PI 
   static member sign n = let k = n % 2 in 1 - double k
   //sum
   member h.sum = function
          | Infinity -> let s = HypHarmonic.sign (m+1) |> float in
                        let p = pown HypHarmonic.pi2 pow in
                        let f = binom.fact pow |> float in
                        s*p*(bernoulli pow) / f |> half
          | Integer i -> [float i..(-1.0)..1.0] |> List.sumBy h.a
   //remainder of a sum
   member h.rest n = h.sum Infinity - h.sum n

//calculate exact distribution for first "apr" terms
//and approximate distribution for the remainder
let distribApprox pow apr x =
  let k = (float (pown 2.0 apr)) in
  let n = Integer apr in                    //number of terms to sum
  let h = HypHarmonic pow in                 //series of even power pow S(n) = 1/1^pow + 1/2^pow + .. +1/n^pow
  let M (h:HypHarmonic) = h.rest n |> half   //remainder of the h's series. It's expectation  Eh
  let m = M h
  let s = M (square h) |> half |> sqrt      //remainder of the h's series. It's stdev  sqrt(Vh)

  //Normal distribution
  let rec erfc x =
      if x < 0.0 then 2.0 - erfc (-x)
      elif x<0.5 then 1.0 - erf x
      elif x>=10.0 then 0.0
      else 
         let P =x*(x*(x*(x*(x*(x*(x*(0.5641877825507397413087057563)
                 + 9.675807882987265400604202961)
                  + 77.08161730368428609781633646)
                   + 368.5196154710010637133875746)
                    + 1143.262070703886173606073338)
                     + 2320.439590251635247384768711)
                      + 2898.0293292167655611275846)
                       + 1826.3348842295112592168999                      
         let Q = x*(x*(x*(x*(x*(x*(x*(1.0 + 17.14980943627607849376131193)
                 + 137.1255960500622202878443578)
                  + 661.7361207107653469211984771)
                   + 2094.384367789539593790281779)
                    + 4429.612803883682726711528526)
                      + 6089.5424232724435504633068)
                       + 4958.82756472114071495438422)
                        + 1826.3348842295112595576438       
         exp (-(square x))*P/Q
  and erf x =
      let z = abs x
      let S = x |> sign |> float
      if z >= 10.0 then S
      else
         if z<0.5 then
          let Xsq = square x 
          let P =Xsq*(Xsq*(Xsq*(Xsq*(Xsq*(Xsq*0.007547728033418631287834
                    + 0.288805137207594084924010)
                     + 14.3383842191748205576712)
                      + 38.0140318123903008244444)
                       + 3017.82788536507577809226)
                        + 7404.07142710151470082064)
                         + 80437.3630960840172832162
          let Q = Xsq*(Xsq*(Xsq*(Xsq*(Xsq + 38.0190713951939403753468)
                     + 658.070155459240506326937)
                      + 6379.60017324428279487120)
                       + 34216.5257924628539769006)
                        + 80437.3630960840172826266
          S*1.1283791670955125738961589031*z*P/Q
         else
            S*(1.0-erfc z) 
            
  let pnormStd t = 
      let sqrt2 = 1.41421356237309504880
      half (erf (t / sqrt2)+1.0) 
  
  let normalize m s x = (x - m) / s in
  //cumulative normal distribution N(m,s)
  let pnorm x m s = x |> normalize m s |> pnormStd 
  
  //builds mediate distribution table for sum S(i-1) + a(i)
  let combine L i =
    let stat = L 
                |> List.fold (fun acc an -> 
                                  let b = h.a (float i) + an in                                
                                  match b with 
                                    | b when b + h.rest (Integer i) < x -> acc        //never reach x -> ignore
                                    | b when b < x -> b::fst acc, snd acc             //lt x -> append to table
                                    | _ -> fst acc, 1.0/float(pown 2.0 i) + snd acc)  //gt x -> add to prob estim
                             (L,0.0)
    //exact table, crude prob estimation
    fst stat, snd stat
   
  //build comlete distribution table for sum S(apr)
  let stat = [1..apr] 
                |> List.fold 
                        (fun acc i -> let c = combine (fst acc) i
                                      (fst c), snd c + snd acc
                        ) ([0.0],0.0)
  
  //compute sum of exact and normal distributions
  fst stat
    |> List.sumBy (fun v -> (1. - pnorm x (m+v) s))  //refine
    |> (*) (inverse k)
    |> (+) (snd stat) 

//comulative distribution for Sum(1,inf, d(i)/i^pow), 
//where P(d(i)=0) = P(d(i)=1) = 0.5 
//tol - converge criteria 10^(-tol) 
let distrib pow tol x = 
   let rec loop cur prev n =
      if abs (prev-cur)>pown 0.1 tol then loop (distribApprox pow n x) cur (n+1)
      else cur
   1. - loop 0.0 1.0 1         

//write to a file   
let write L =
     System.IO.File.WriteAllLines(@"c:\temp\a.csv" ,L)
     printfn "OK"

//write out cumulative distribution points to a file 
[|0.0..0.001..1.7|] |> Array.map (fun x -> sprintf "%f, %f" x (distrib 2 8 x)) |> write
val square : x:'a -> 'b (requires member ( * ))

Full name: Script.square
val x : 'a (requires member ( * ))
Multiple items
val double : x:'a -> 'b (requires member ( + ))

Full name: Script.double

--------------------
type double = System.Double

Full name: Microsoft.FSharp.Core.double
val x : 'a (requires member ( + ))
val half : x:float -> float

Full name: Script.half
val x : float
val inverse : x:float -> float

Full name: Script.inverse
val sign : x:float -> float

Full name: Script.sign
val abs : value:'T -> 'T (requires member Abs)

Full name: Microsoft.FSharp.Core.Operators.abs
type binom =
  static member comb : n:int -> m:int -> int
  static member fact : n:int -> int

Full name: Script.binom
static member binom.fact : n:int -> int

Full name: Script.binom.fact
val n : int
val loop : (int -> int -> int)
val acc : int
static member binom.comb : n:int -> m:int -> int

Full name: Script.binom.comb
val m : int
static member binom.fact : n:int -> int
val bernoulli : _arg1:int -> float

Full name: Script.bernoulli
Multiple items
module List

from Microsoft.FSharp.Collections

--------------------
type List<'T> =
  | ( [] )
  | ( :: ) of Head: 'T * Tail: 'T list
  interface IEnumerable
  interface IEnumerable<'T>
  member GetSlice : startIndex:int option * endIndex:int option -> 'T list
  member Head : 'T
  member IsEmpty : bool
  member Item : index:int -> 'T with get
  member Length : int
  member Tail : 'T list
  static member Cons : head:'T * tail:'T list -> 'T list
  static member Empty : 'T list

Full name: Microsoft.FSharp.Collections.List<_>
val map : mapping:('T -> 'U) -> list:'T list -> 'U list

Full name: Microsoft.FSharp.Collections.List.map
val k : int
static member binom.comb : n:int -> m:int -> int
Multiple items
val float : value:'T -> float (requires member op_Explicit)

Full name: Microsoft.FSharp.Core.Operators.float

--------------------
type float = System.Double

Full name: Microsoft.FSharp.Core.float

--------------------
type float<'Measure> = float

Full name: Microsoft.FSharp.Core.float<_>
val sumBy : projection:('T -> 'U) -> list:'T list -> 'U (requires member ( + ) and member get_Zero)

Full name: Microsoft.FSharp.Collections.List.sumBy
val x : float * float
val fst : tuple:('T1 * 'T2) -> 'T1

Full name: Microsoft.FSharp.Core.Operators.fst
val snd : tuple:('T1 * 'T2) -> 'T2

Full name: Microsoft.FSharp.Core.Operators.snd
type UpperLimit =
  | Infinity
  | Integer of int

Full name: Script.UpperLimit
union case UpperLimit.Infinity: UpperLimit
union case UpperLimit.Integer: int -> UpperLimit
Multiple items
val int : value:'T -> int (requires member op_Explicit)

Full name: Microsoft.FSharp.Core.Operators.int

--------------------
type int = int32

Full name: Microsoft.FSharp.Core.int

--------------------
type int<'Measure> = int

Full name: Microsoft.FSharp.Core.int<_>
Multiple items
type HypHarmonic =
  new : pow:int -> HypHarmonic
  member a : i:float -> float
  member pow : int
  member sum : (UpperLimit -> float)
  member rest : n:UpperLimit -> float
  static member pi2 : float
  static member ( * ) : a:HypHarmonic * b:HypHarmonic -> HypHarmonic
  static member sign : n:int -> int

Full name: Script.HypHarmonic

--------------------
new : pow:int -> HypHarmonic
val pow : int
val repl : (float -> float)
val loop : (('a -> 'a) -> int -> 'a -> 'a)
val f : ('a -> 'a)
val h : HypHarmonic
member HypHarmonic.pow : int

Full name: Script.HypHarmonic.pow
member HypHarmonic.a : i:float -> float

Full name: Script.HypHarmonic.a
val i : float
val a : HypHarmonic
val b : HypHarmonic
property HypHarmonic.pow: int
static member HypHarmonic.pi2 : float

Full name: Script.HypHarmonic.pi2
namespace System
type Math =
  static val PI : float
  static val E : float
  static member Abs : value:sbyte -> sbyte + 6 overloads
  static member Acos : d:float -> float
  static member Asin : d:float -> float
  static member Atan : d:float -> float
  static member Atan2 : y:float * x:float -> float
  static member BigMul : a:int * b:int -> int64
  static member Ceiling : d:decimal -> decimal + 1 overload
  static member Cos : d:float -> float
  ...

Full name: System.Math
field System.Math.PI = 3.14159265359
static member HypHarmonic.sign : n:int -> int

Full name: Script.HypHarmonic.sign
member HypHarmonic.sum : (UpperLimit -> float)

Full name: Script.HypHarmonic.sum
val s : float
static member HypHarmonic.sign : n:int -> int
val p : float
val pown : x:'T -> n:int -> 'T (requires member get_One and member ( * ) and member ( / ))

Full name: Microsoft.FSharp.Core.Operators.pown
property HypHarmonic.pi2: float
val f : float
val i : int
member HypHarmonic.a : i:float -> float
member HypHarmonic.rest : n:UpperLimit -> float

Full name: Script.HypHarmonic.rest
val n : UpperLimit
property HypHarmonic.sum: UpperLimit -> float
val distribApprox : pow:int -> apr:int -> x:float -> float

Full name: Script.distribApprox
val apr : int
val k : float
val M : (HypHarmonic -> float)
member HypHarmonic.rest : n:UpperLimit -> float
val m : float
val sqrt : value:'T -> 'U (requires member Sqrt)

Full name: Microsoft.FSharp.Core.Operators.sqrt
val erfc : (float -> float)
val erf : (float -> float)
val P : float
val Q : float
val exp : value:'T -> 'T (requires member Exp)

Full name: Microsoft.FSharp.Core.Operators.exp
val z : float
val S : float
val Xsq : float
val pnormStd : (float -> float)
val t : float
val sqrt2 : float
val normalize : (float -> float -> float -> float)
val pnorm : (float -> float -> float -> float)
val combine : (float list -> int -> float list * float)
val L : float list
val stat : float list * float
val fold : folder:('State -> 'T -> 'State) -> state:'State -> list:'T list -> 'State

Full name: Microsoft.FSharp.Collections.List.fold
val acc : float list * float
val an : float
val b : float
val c : float list * float
val v : float
val distrib : pow:int -> tol:int -> x:float -> float

Full name: Script.distrib
val tol : int
val loop : (float -> float -> int -> float)
val cur : float
val prev : float
val write : L:string [] -> unit

Full name: Script.write
val L : string []
namespace System.IO
type File =
  static member AppendAllLines : path:string * contents:IEnumerable<string> -> unit + 1 overload
  static member AppendAllText : path:string * contents:string -> unit + 1 overload
  static member AppendText : path:string -> StreamWriter
  static member Copy : sourceFileName:string * destFileName:string -> unit + 1 overload
  static member Create : path:string -> FileStream + 3 overloads
  static member CreateText : path:string -> StreamWriter
  static member Decrypt : path:string -> unit
  static member Delete : path:string -> unit
  static member Encrypt : path:string -> unit
  static member Exists : path:string -> bool
  ...

Full name: System.IO.File
System.IO.File.WriteAllLines(path: string, contents: System.Collections.Generic.IEnumerable<string>) : unit
System.IO.File.WriteAllLines(path: string, contents: string []) : unit
System.IO.File.WriteAllLines(path: string, contents: System.Collections.Generic.IEnumerable<string>, encoding: System.Text.Encoding) : unit
System.IO.File.WriteAllLines(path: string, contents: string [], encoding: System.Text.Encoding) : unit
val printfn : format:Printf.TextWriterFormat<'T> -> 'T

Full name: Microsoft.FSharp.Core.ExtraTopLevelOperators.printfn
module Array

from Microsoft.FSharp.Collections
val map : mapping:('T -> 'U) -> array:'T [] -> 'U []

Full name: Microsoft.FSharp.Collections.Array.map
val sprintf : format:Printf.StringFormat<'T> -> 'T

Full name: Microsoft.FSharp.Core.ExtraTopLevelOperators.sprintf
Next Version Raw view Test code New version

More information

Link:http://fssnip.net/7Xw
Posted:4 years ago
Author:Pavel Tatarintsev
Tags: combinatorics , normal distribution , bernoulli numbers , hyperharmonic series