(*Edgeworth expansion for Stirling numbers of the first kind
We use the expansion and the notation from the paper 
"General Edgeworth expansions with applications to profiles of random trees"
 by Zakhar Kabluchko, Alexander Marynych, Henning Sulzbach
available at http://arxiv.org/abs/1606.03920
*)

(*Parameters: 
phi(beta) = Exp[beta]-1,
  LogW[beta] = LogGamma[theta]-LogGamma[theta*Exp[beta]],
  kappa[[j]] = j-th derivative of phi(beta)=Exp[beta]-1  at beta=0,
  chi[[j]] = j-th derivative of LogW at beta=0,
  sigma[[j]] = 1 is ignored everywhere.
 *)
kappa=Table[1,{n,1,5}]
LogW[beta_] = LogGamma[theta]-LogGamma[theta*Exp[beta]]
chi=List[D[LogW[x],x]/.x->0,D[D[LogW[x],x],x]/.x->0, D[D[D[LogW[x],x],x],x]/.x->0]

(*Differential operators D_1,D_2,D_3. Here D stays for sigma^(-1) d/dx *)

D1= kappa[[3]]/6*D^3 + chi[[1]]*D^1
D2= kappa[[4]]/12*D^4+ chi[[2]]*D^2
D3= kappa[[5]]/20*D^5 + chi[[3]]*D^3

(*
Bell polynomials are 
B_0=1,
 B_1(z_1) = z_1,
 B_2(z_1,z_2) = z_1^2 + z_2,
 B_3(z_1,z_2,z_3) = z_1^3 + 3z_1z_2 + z_3
*)
(*
Probabilist Hermite polynomials are given by
He_p(x) = 2^(-p/2)*HermiteH[p, x/Sqrt[2]]
To check this, use: Table[Simplify[2^(-p/2)*HermiteH[p, x/Sqrt[2]]], {p,0,6}]
*)

(*Terms in the Edgeworth expansion of the Stirling numbers are denoted by G_0(x), G_1(x),...*)
H0=1;

H1[x_]=Simplify[ (Expand[D1]/. D^p_ -> 2^(-p/2)*HermiteH[p,x/Sqrt[2]]) /. D -> 2^(-1/2)*HermiteH[1,x/Sqrt[2]]]

A11=Simplify[Coefficient[H1[x], x]]
A12=Simplify[Coefficient[H1[x], x^3]]

H2[x_]=Simplify[1/(2!)*Expand[D1^2 + D2]/. D^p_ -> 2^(-p/2)*HermiteH[p,x/Sqrt[2]]]

A21=Simplify[H2[0]]
A22= Simplify[Coefficient[H2[x], x^2]]

H3[x_]=Simplify[1/(3!)*Expand[D1^3 + 3 D1*D2+D3]/. D^p_ -> 2^(-p/2)*HermiteH[p,x/Sqrt[2]]]

A31= Simplify[Coefficient[H3[x],x]]

(*
The first three terms in the expansion for x = a/Sqrt[w]. 
The term H4/w^2 = constant/w^2 +o(1/w^2) is ignored because the constant does not depend on a.  
Error: o(1/w^2)
*)
Series[Exp[-t^2/2], {t,0,4}]

A= (1+H1[a/Sqrt[w]]/Sqrt[w] + H2[a/Sqrt[w]]/w +H3[a/Sqrt[w]]/w^(3/2)) * (1-a^2/(2w) + 1/8 *a^4/w^2);
Collect[Coefficient[A,1/w],a]
Collect[Coefficient[A,1/w^(1/2)],a]
Collect[Coefficient[A,1/w^(3/2)],a]
Collect[Coefficient[A,1/w^2],a]

(*Here we insert a = gam + chi[[1]]-1/2, where gam is the new variable*)Poly[gam_]= Simplify[Coefficient[A,1/w^2] /. a->gam+chi[[1]]-1/2]

(*
Now we can compute the quantity denoted by s^*(theta) in the paper "Asymptotic expansions for the Ewens distribution and the Stirling numbers of the first kind"
*)
Simplify[Poly[1/2] -Poly[-1/2]]

(*The resulting expression is positive*)
Plot[Poly[1/2] -Poly[-1/2], {theta,0,1}]