May be, that 'tryhard' has problems with uninitialized variables?
I do not know, but simplify/size is quite good - I included that.

You do not always need 'tryhard': if it is OK without than it is
OK. And do not forget, that finally your (Fortran) compiler will
also do an optimization, so it will be altered in any case.

Also the dimensions of arrays FORTR and stampa do not fit for
assigning, just correct it.

The uploaded sheet contains the modifications, searching for my
initials "AVt" shows some comments.

For the question "codegen or CodeGeneration?": if Maple kills all
the old stuff (at least from the help pages), then it really should
provide equivalents - and they dont (makeproc). What sloppiness.

Hope it helps ...

Download 102_codegen_Fortran_Kram.mws
View file details

More or less - if one wants it numerical (since there seems to be
no easy chance for an analytic solution ... but it is only a model,
to be justified, so perhaps a modification would do as well ...),
then more has to invested (due to numerical cancellation problems
through extreme magnitudes).

Here is a sketch: one can do the integration wrt z (Maple seems to
hang up in simplifications, due to concrete figures), contributing
a multiplicative constant (in terms of Gamma).

The result is like y^r*(-y+1)^(-s)*F*E, where E=exp(quadric in x)
and F is the rest (powers + log).

Pulling out the constant of the quadric it can be seen to be affine
equivalent to sum( x_i ^2 ) + some constant (i.e. a sphere).

Thus (up to constants) we have y^r*(-y+1)^(-s)*F*exp(sum( x_i ^2 ))
by a multi-dimensional change of coordinates in x (note that one
need not to care for domains, since affines map IR^n onto itself
and the needed det(Jacobian) is a constant [y is not affected]).

The last term is a product of normal distributions, hence centered
at zero - which gives a guess, where to integrate: the normal pdf's
should dominate the other terms.

Now carefully collecting all constants *may* give a chance, even if
it seems the integrand is quite small ( << 1).

But I do not want to carry out that plan. And may be some further
simplifications are possible (looking at 'F' and y and log(sum exp)).

However my objection is, that if it is already that difficult to
find the normalizing constant (and that might be somewhat incorrect),
then a subsequent use has a chance to be numerical difficult as well.


PS: first though the poster is actually 'Mike' ... that would be
typical for him ...

1st: avoid floats like C=1.000 if you have no real reason, try exact values.

2nd: try do use functions, it gives better structure

3rd: if having doubts and using test values, then try simple ones or just one
(which can be typed in without running a script).

4th: avoid 'displayprecision' if there is no need for it, try evalf[3]().

Hm ... and using Array (upper case) is more concurrent ...


Now again:

  a := 1: C := 1: L := 4: CS := C*L: # avoid floats ...
  
  # your formula:
  Product(
    GAMMA(1+CS)*GAMMA(1+CS+(L-n+1)*r)*Hypergeom([1+CS, 1], [1+C], r)/
      (GAMMA(r+1+CS)*GAMMA(1+CS+(L-n)*r)), n = 1 .. L-1)*
    Hypergeom([1+CS, 1], [1+C], r);  
  #simplify(%); # not really needed, can be done later ...
  value(%);
  convert(%,StandardFunctions); 
  factor(%);
  f:=unapply(%,r);

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


Your values:

  theR:=Array(1 .. 10):
  
  ind := 1:
  for rangedb from -5 by 5 to 40 do 
    y := 10^((1/10)*rangedb); 
    R := C/y; 
    r := R/(a+L*R); #r:='r':
    theR[ind] := r; 
    ind := ind+1 
  end do:
  r:='r':

Do not forget to reset your variables ... I would do it another way.

  theR:=evala(theR);

This gives the exact values for your various r. For all of them you
claim that different results are to be given. So take one, the 2nd
is quite simple, it is r = 1/5.

  f(1/5);        # strange, what Maple does here? No, due to Gamma.
  simplify(%);
  evalf[200](%): # use high precision in case of doubts
  evalf(%);      # cut down do working precision

                             11.07390057

Thx for the paper, however I do not see your formula and the notations
in the article are quite technical, and I do not want to step into the
details there. Well, I do not even see, which coefficients and inputs
you have to use F_A in eq (1) of the paper.

However I am sure your code is not correct, i.e. Product()* 2F1 is not
what you want.


For the convergence problem (approaching abs = 1): this already occures
if working with hypergeometric series, thus Maples does not use it in
a direct, simple way.

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.