To get the above expression I first computed your original
integrand 'PostNumG1' using 100 Digits, converted floats to
Rationals. It turns out, that an integration over x4 ... x8
can be done symbolically and leads to 'the' K above.

To speed up computation for your integrand one can use
obvious assumptions for the 2 procedures, which use some
'fsolve' (positivity).

However I do not see a way for much more symbolic integration.

Thx, I also meant in which field all that should be applied ...

I think the approach will not work in usual double precision
and thus method = _cuhre will not work.

After some brute roundings you can roughly approximate your
integral byit using

  2*10^12 * exp(-20000) * Int(smallerK, x,y,z), where

smallerK:= (-z+1)^(4/5)/(-y+1)^(2/25)*(z-1)^10*y^20*z^40*
  (y*exp(x1)+y*exp(x2)+exp(x3)+exp(x4))^900 *
  exp(-500*x1^2+2*x1*x2+900*x1*x3+2*x1*x4-4/5*x2^2+
       2*x2*x3+3/10*x2*x4-500*x3^2+2*x3*x4-4/5*x4^2+
       1000*x1-30*x2+1000*x3-30*x4)*
  exp(-50*ln(y*exp(x1)+y*exp(x2)+exp(x3)+exp(x4))^2);

  x = (x1,x2,x3,x4) in IR^4, y and z are in ]0 ... 1[.

Due to your description that integral should be non-zero,
at least analytical this is true.

Shooting some values for that integrand (or your original)
shows, that it quickly leads to numerical overflow/underflow,
even for x=0, y=z=1/2.

I do not see in which region(s) it would give reasonable and
usable values, but it seems to be 'small' compared to IR^4.

This would mean, that it would become strongly peaked, if
one uses a transform IR^4 to ]0 ...1[^4.

So in any case you need to locate a region, where to work
and need reasonings why off that the integration error can
be neglected. Good luck ...

Running the inner loop with uppercase Hypergeom, an unsigned r and CS=L*C one gets
a result, which lets one guess the original formula used here (by varying L):

hypergeom([1, 1+L*C],[1+C],r)^L *
  GAMMA(1+(C+r)*L)/(GAMMA(r+1+L*C)^(L-1))/GAMMA((C+r)*L+1-3*r)*GAMMA(1+L*C)^(L-1);

Substituting postive integers for L and C this turns out to be a rational
function in r, up to factors involving the Gamma function:

  subs(L=4, C = 1, %);
  convert(%,StandardFunctions):
  factor(%):
  f:=unapply(%,r); 

                                            4   2           4
                   54 GAMMA(5 + 4 r) (r - 2)  (r  - 2 r + 2)
         f := r -> ------------------------------------------
                                   16             4
                            (r - 1)   GAMMA(r + 5)


The outer loop only defines test values, which may be sampled into an Array and
applies f to that. Or one does it follows:

  data:=[a = 1, C = 1, L = 4];

  for rangedb from -5 by 5 to 40 do 
    y := 10^((1/10)*rangedb); 
    R := C/y; 
    r := R/(a+L*R);
    eval(f(r),data);
  print('f'(r)=evalf[20](%)); #print('f'(eval(r,data))=evalf[20](%));
  end do:

Whether those are results coinciding with the article now is on you, since that
paper is not public available :-)

For your series problem: using series for hypergeometrics (if even allowed to
be used) becomes quite a mess for inputs x where abs(x) > 2/3 and Maple uses
other ways. However usually you do not need to care how Maple does it, just
call the function.

What is the magnitude of your expected result? And again: Where
does that problem come from, for what should the result be used?

Retry it using a classical sheet (file extension *.wms) and upload it
(using the 'light, green arrow' in the editor for replying mails here).

If your Sigma is numerical, then the MNormDist becomes a (real) constant.

You have PostNumG1 = A * E (see above). But these terms to not contain
the all variables, which you use for integration (MNormDist = constant),
except I understood your sheet wrong and made some error.

But this may be not what you want: assuming that also sigma[G] is some
constant (not 0) then the inner integral is A * Int(E) over IR^4, and
that Int(E) is a function f(y) of y only.

Thus you want Int( A * f(y) ) using x[5], ..., x[8], y, z. 

But already the inner integral here is infinity, since the A*f(y) is
some constant w.r.t. the x[5], so over the full Reals integral=infinity.

For me it seems you have set up some wrong formula, may be through some
typos or handling problems.

So it is more than a 'numerical issue' ...

Edited to add:

And as already said: your variables are LOWER case x, while you use UPPER case X
in that evalf ...

The term E will give complex results: using t=x[1] and s= sigma[G] it writes like
exp(-1/8*(2*ln(a*t+b)+s^2)^2/s^2) for a and b real constants. Now integrating over
the Reals w.r.t. the t should result in complex numbers due to log on negatives.

Moreover one sees, that this is not a (multivariate) normal.


Where does that problem come from, for what should the result be used?

I am not sure, whether your task actually evaluates to a number (so one
has to choose a numeric integration method), since you post it in symbols
to give us an overview.

The posted coded will hang up at 

  a3:=1/Sigma;        # since it means to invert an 8 x 8 symbolic matrix

  MNormDist := 1/((2*Pi)^4*sqrt(Determinant(Sigma)))*a4: # similar reason

Omitting 'evalf' shows the integrand as A*E using

  tmp:=PostNumG1; #tmp:=op(1,DenomG1); 
  E:=select(has,%,exp); lprint(%);

    exp(-1/2*(ln(20000)-ln(y*(x[1]+x[2])+x[3]+x[4])-1/2*sigma[G]^2)^2/sigma[G]^2)

  A:=tmp/E; lprint(%);

    1/40000*MNormDist*y^(-1+AlphaA)*(1-y)^(BetaA-1)/Beta(AlphaA,BetaA)*z^(-1+AlphaB)*
      (1-z)^(BetaB-1)/Beta(AlphaB,BetaB)/sigma[G]*2^(1/2)/Pi^(1/2)

Now 

  PostNumG1; indets(%,atomic);

shows you need to integrate over x[i], not X[i] (lower case vs uppercase) and it
does *not* contain x[5] ... x[8] (except they are hidden in MNormDist or similar).

For this also note, that A does not depend on x[1] ... x[4] (except MNormDist?).


So formally it is a bit odd ... may be due to the intended simplification to
give a convenient overview.


For the omitted steps (a3 and MNormDist) I would at least high(er) precision and
cut down later.


For the actual numerical integration ... well ... äh ...

It is more likely you have typos in you code or translated the formula
not in a correct way - however can not say something on that (do not
work with that functions and do not have the article).

But looking into your example: the outer loop just produces different
inputs r ( 0 < r < 1/4) for the inner loop.

The inner loop leads to products of Gamma and powers of 2F1. Even if
some cancellation error is possible (may be healed by high precision)
I doubt that Maple makes an error in that case.


PS: it is better to work with the classical interface for such kind of
problems, there are less traps for false input (say multiplication is
missing or similar).

Edited: you inner (with the given parameters) loop produces
54*(r-2)^4*(r^2-2*r+2)^4*GAMMA(5+4*r)/(r-1)^16/GAMMA(r+5)^4

not a structural suggestion (do not have a good image or idea for your problem):

how about using 'evalhf' or even 'Compile' (as it seems you want it of that exactness) ?

You are supposed to

1. Use Maple syntax (as nobody likes to type in your problem)

2. Ask a question ("any advice?" might be answered by "sure, have a coffee!")

3. It does not hurt to speak in complete sentences

And in most places some 'magic word' is not a bad idea ...

So please give your problem a new chance, yes?

Not that I am going to prove it, but if the inverse is dense,
then I would guess it may be a quite complicated expression.

Here is a way to check it on test values:

  restart; 
  interface(rtablesize=14); with(LinearAlgebra); Digits:=14;

  m:=RowDimension(M);
  indets(M): L:=convert(%,list); k:=nops(%);

  # the Matrix is of dim=14, you have 27 variables
  # set them to some values

  R:=RandomTools[Generate](list( integer(range=-20..20), k));
  data:=[seq( L[i]=R[i], i=1..k)];

  N:=eval(M,data); 
  Q:=MatrixInverse(N);

Then Q is returned almost immediately - and it is dense.

And(And(abs(a)<infinity,a<>0),And(abs(b)<infinity,b<>0));
subs(And=`and`,%);

  and(and(| a | < infinity, a <> 0), and(| b | < infinity, b <> 0))

piecewise(%,1,0);

                     { 1        a <> 0 and b <> 0
                     {
                     { 0            otherwise

Hm ... it is not so much fun to check codes - and yours is extremely ugly (sorry to say so, it is not meant personal offending).

You permanently overwrite variables, no line breaks, no idents ...

My guessing is: you should reset the variables for the loops (or put the inner body into a procedure, which would force you to be more clean).

The error message gives you, where the first error occurs: when you use 'lhs'

Hope you can accept, that I do not want to step in details through your code - may be, the suggestions above give you some idea, how to build it more clean - otherwise it is hard to check.

If the question is as I suppose (but that guy seems to be too lazy
to answer the implicite question), then I think there is not extremum
except at the boundaries.

E0:=Int(q*exp((1-b)^2* z^2/lambda )/((1+b^2*q^2+b^2*z^2)*(q^2+z^2)),
  z=-infinity..infinity);

E:=eval(E0, lambda=-alpha^2); # since 0 < lambda ===> infinite
E1:=PDETools[dchange](z = x/(1-b)*alpha, %, [x]) assuming 0 < alpha:

That depends on sign( b-1 ), but changes only the sign (use 'Flip').

S:=[b = 1/(s^(1/2)*beta), q = p^(1/2)*beta, alpha = beta-1/(s^(1/2))];

  E1 assuming b-1 < 0;
  eval(%, S);
  PS:=simplify(%) ;

                       infinity
                      /                    2   1/2
                     |             s exp(-x ) p
              PS :=  |          --------------------- dx
                     |                2    2
                    /           (p + x ) (x  + p + s)
                      -infinity

  plot3d(PS, s=0..2,p=20..1, axes=boxed);

This plot suggests, that the value is decreasing with increasing p
(but I am too lazy too for trying a proof), which should give that 
there is no finite minimum.

I think he wants to minimize an expectation value w.r.t. that q
and that the integral (up to factors) can be written as 

  Int( exp(-x^2) * f(x,q), x = -infinity .. infinity), 
  f(x,q)=q/polynom(x,q)

where that polynomial in x has 4 distinct complex roots.

Maple does not give the integral in closed form (even after using
(real) partial fractions), i.e. the special function(s) involved
are possibly not of known form. And even if it would, it is not
clear that if would lead to some usefull result for q_min.

And 0 = diff(theIntegral,q) looks not that nice to solve ...

May be VariationalCalculus can be used. But I guess, one needs
sound Math, not only Maple, for a usefull description of the
desired aim.

I always forget the 'distributed' option for collect ...

One can also use coeftayl(p,c0=0,1) or coeff(p,c0,1)