1)Use HGM and NIntegrate to compute the 4-fold  integral  of the expection of an Euler characteristic number

<< RISC`HolonomicFunctions` (* load Christoph's package *)

HolonomicFunctions Package version 1.7.3 (21-Mar-2017)
written by Christoph Koutschan
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

ann3 = << ec1_ann3.m; (* input the annihilator of the inner triple integral of Example 1*)

uts3 = UnderTheStaircase[{ann3}] (* the parametric derivatives under the stair of ann3 *)

stheta=(2 s)/(1+s^2); (* make a universal transformation for Sin[theta] *)

ctheta=(1-s^2)/(1+s^2); (* make a universal transformation for Cos[theta] *)

sphi=(2 t)/(1+t^2); (* make a universal transformation for Sin[phi] *)

cphi=(1-t^2)/(1+t^2);(* make a universal transformation for Cos[phi] *)

R=s1 (b sphi stheta+d cphi ctheta-m11)^2+s2 (-b sphi ctheta+d cphi stheta-m21)^2+s2 (b cphi ctheta+d sphi stheta-m22)^2+s1 (-b cphi stheta+d sphi ctheta)^2; (* input R in the integrand *)

f=2/(1+t^2 )2/(1+s^2) (d^2-b^2) (s1 s2)/(2 Pi)^2 Exp[-1/2 R]; (*input the integrand of Example 1*)

f1=Simplify[f/.{s1->2,s2->1,m11->1,m21->-1,m22->1}] ;(* specify the parameters as that in Example 1*)

Timing[inits = NIntegrate[Together[ApplyOreOperator[uts3,1/2 f1] /. {d -> 1/2}], {b, -Infinity, Infinity}, {s, -Infinity, Infinity}, {t, -Infinity, Infinity},MaxRecursion->20,Method->{GlobalAdaptive,MaxErrorIncreases->10000}];] (* compute the initial value for the Pfaffian system of ann3 at d = 1/2. Since ann3 has a singularity at d = 0, it is not suitable to choose the initial point at d = 0 *)

eqnsd = ApplyOreOperator[ann3, y[d]] == 0; (* set up the differential equation for the inner triple integral *)

inits1 = Table[D[y[d], {d, i}] == inits[[i +1 ]], {i, 0, 9}] /. {d -> 1/2} (* set up the initial values for eqnsd at d = 1/2 *)

s1 = NDSolve[Join[{eqnsd},inits1], y, {d, 1/2, 30}]; (* numerical solving of the diffential equation eqnsd from 1/2 to 30 *)

Plot[Evaluate[y[d]/.s1],{d,1/2,30},PlotRange->All] (* y[d] is very small when d ≥ 5*)

HGM1 = Table[NIntegrate[y[d] /. s1, {d, i, 30}][[1]], {i, 1, 6}] (* Using numerical integration to compute the integral in Example1 for x from 1 to 6. The results are as accurate as that in Prof. Kuriki's mathematica notebook*)

2) Monte Carlo study for the Expectation of Euler charateristic number  of the expection of an Euler characteristic number

EulerChar[A_,x_]:=Module[{b,d,a,c},{b,d}=SingularValueList[A];(*b and d are singular value of A*)If[TrueQ[b≥x],a=Sign[b^2-d^2],a=0];
If[TrueQ[d≥x],c=Sign[d^2-b^2],c=0];
Return[a+c];]; (* compute the Euler charateristic number of a random matrix A w.r.t. x*)

Clear[s1]

sub={s1->2,s2->1,m11->1,m21->-1,m22->1}

S=DiagonalMatrix[{1/Sqrt[s1],1/Sqrt[s2]}]/.sub

M={{m11,0},{m21,m22}}/.sub

n =10000000 (* the number of iterations *)

ra=Table[S.{RandomVariate[NormalDistribution[],2,WorkingPrecision->10],RandomVariate[NormalDistribution[],2,WorkingPrecision->10]}+M,{i,1,n}]; (* generate n random matrix A*)

Timing[mc = Table[ N[Sum[EulerChar[ra[[i]],x],{i,1,n}]/n], {x, 1, 6}];] (*the results of expectation of Euler characteristic number by Monte Carlo study *)

mc (* The results are as accurate as that in Prof. Kuriki's mathematica notebook *)