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 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