Your sheet has some typos, and I guess you assume 0<y, 0<e, but j(e,y) is unclear.
In your sheet you have e^something, but I guess you mean exp(something), right?

Then I think you want the following

dummy := sqrt((e*y+1)/(y*e));
dummy1 := sqrt(1/(e*y+1));
Q := RandomVariable(Normal(0, dummy));
R := RandomVariable(Normal(y*e*z/(e*y+1), dummy1));

'Int(Int(-exp(-a*j(e, y)*r)*PDF(R, r), r = -infinity .. infinity)*PDF(Q, z), 
  z= -infinity .. infinity)';

has a solution -exp(1/2*j(e,y)^2*a^2), so its maximum is the minimum of j(e,y).

The 2 other integrals also have explicit solutions of type pdfN*cdfN.

Hence you do not need any Gauss-Hermite but have to solve your task from there.

PS: do not use lower case "int", use Uper Case "Int" and the command "value".

PPS: for beginners I would suggest to use the 'Classical Interface' (like my
appended sheet). It does not look that nice - but quite often it is more safe
to be used and clearly shows what one has typed in. A matter of taste, yes.

MP_147868-GaussHermi.mws

Your sheet has some typos, and I guess you assume 0<y, 0<e, but j(e,y) is unclear.
In your sheet you have e^something, but I guess you mean exp(something), right?

Then I think you want the following

dummy := sqrt((e*y+1)/(y*e));
dummy1 := sqrt(1/(e*y+1));
Q := RandomVariable(Normal(0, dummy));
R := RandomVariable(Normal(y*e*z/(e*y+1), dummy1));

'Int(Int(-exp(-a*j(e, y)*r)*PDF(R, r), r = -infinity .. infinity)*PDF(Q, z), 
  z= -infinity .. infinity)';

has a solution -exp(1/2*j(e,y)^2*a^2), so its maximum is the minimum of j(e,y).

The 2 other integrals also have explicit solutions of type pdfN*cdfN.

Hence you do not need any Gauss-Hermite but have to solve your task from there.

PS: do not use lower case "int", use Uper Case "Int" and the command "value".

PPS: for beginners I would suggest to use the 'Classical Interface' (like my
appended sheet). It does not look that nice - but quite often it is more safe
to be used and clearly shows what one has typed in. A matter of taste, yes.

MP_147868-GaussHermi.mws

Your 1st (for z <= 0) and 2nd integral (for z < 0) can be given formally as

Int1:= (e, j, y, z) -> -1/2/(2*e*y+2)^(1/2)*(-erf(1/8/(2*e*y+2)^(1/2)*(8*e*y*z-
3*j))+1)*(e*y+1)^(1/2)*2^(1/2)*exp((-48*e*j*y*z+9*j^2)/(128*e*y+128))-1/2*2^(1/2)*
(e*y+1)^(1/2)*(erf(1/(2*e*y+2)^(1/2)*y*e*z)+1)/(2*e*y+2)^(1/2);

Int2:= (e, j, y, z) -> -1/2*2^(1/2)*(e*y+1)^(1/2)*exp((-48*e*j*y*z+9*j^2)/(128*e*y
+128))*(erf(1/8/(2*e*y+2)^(1/2)*(8*e*y*z-3*j))+1)/(2*e*y+2)^(1/2)-1/2*2^(1/2)*(e*y+1)
^(1/2)*(-erf(1/(2*e*y+2)^(1/2)*y*e*z)+1)/(2*e*y+2)^(1/2);

Using test values tst1:=[y=1, e=0.1e-1, z=h[19]] and  tst2:=[y=1, e=0.1e-1, z=h[1]];
it seems that these do not have a maximum for finite j in IR.

Generally you have y::posint, 0 <= e and from that (and 0 <= z or z < 0) one can
compute the limit of both infinities

L1 := 1/2*(-2^(1/2)*(e*y+1)^(1/2)*erf(1/(2*e*y+2)^(1/2)*y*e*z)-2^(1/2)*
(e*y+1)^(1/2))/(2*e*y+2)^(1/2);
L2 := 1/2*(2^(1/2)*(e*y+1)^(1/2)*erf(1/(2*e*y+2)^(1/2)*y*e*z)-2^(1/2)*
(e*y+1)^(1/2))/(2*e*y+2)^(1/2);

So I think that using the Optimizer returns a faked result (due to coarse tolerance).

I think there are only finite minima, these are the values for j such that INT = 0 .

Your 1st (for z <= 0) and 2nd integral (for z < 0) can be given formally as

Int1:= (e, j, y, z) -> -1/2/(2*e*y+2)^(1/2)*(-erf(1/8/(2*e*y+2)^(1/2)*(8*e*y*z-
3*j))+1)*(e*y+1)^(1/2)*2^(1/2)*exp((-48*e*j*y*z+9*j^2)/(128*e*y+128))-1/2*2^(1/2)*
(e*y+1)^(1/2)*(erf(1/(2*e*y+2)^(1/2)*y*e*z)+1)/(2*e*y+2)^(1/2);

Int2:= (e, j, y, z) -> -1/2*2^(1/2)*(e*y+1)^(1/2)*exp((-48*e*j*y*z+9*j^2)/(128*e*y
+128))*(erf(1/8/(2*e*y+2)^(1/2)*(8*e*y*z-3*j))+1)/(2*e*y+2)^(1/2)-1/2*2^(1/2)*(e*y+1)
^(1/2)*(-erf(1/(2*e*y+2)^(1/2)*y*e*z)+1)/(2*e*y+2)^(1/2);

Using test values tst1:=[y=1, e=0.1e-1, z=h[19]] and  tst2:=[y=1, e=0.1e-1, z=h[1]];
it seems that these do not have a maximum for finite j in IR.

Generally you have y::posint, 0 <= e and from that (and 0 <= z or z < 0) one can
compute the limit of both infinities

L1 := 1/2*(-2^(1/2)*(e*y+1)^(1/2)*erf(1/(2*e*y+2)^(1/2)*y*e*z)-2^(1/2)*
(e*y+1)^(1/2))/(2*e*y+2)^(1/2);
L2 := 1/2*(2^(1/2)*(e*y+1)^(1/2)*erf(1/(2*e*y+2)^(1/2)*y*e*z)-2^(1/2)*
(e*y+1)^(1/2))/(2*e*y+2)^(1/2);

So I think that using the Optimizer returns a faked result (due to coarse tolerance).

I think there are only finite minima, these are the values for j such that INT = 0 .

Please avoid any "conditional" and use the formal expressions

Please avoid any "conditional" and use the formal expressions

Hm ... you commented the code.

What I meant: write down the explicit integrals and the task (without using the nodes,
which would Numerics) in terms of Analysis (not Statistics), then what should be 'optimized'

PS: h := evalf(sqrt(2)*ysigma*p) does not make (formal) sense, your p is an Array
but you use h as a Real - a missing indexing ?

Hm ... you commented the code.

What I meant: write down the explicit integrals and the task (without using the nodes,
which would Numerics) in terms of Analysis (not Statistics), then what should be 'optimized'

PS: h := evalf(sqrt(2)*ysigma*p) does not make (formal) sense, your p is an Array
but you use h as a Real - a missing indexing ?

It is pure code and hard to read and to understand. This way I do not see,
why you need Gauss-Hermite.

Could you drop some explanations? For the actual task, if possible, and
the method used

At least both the integrals can be computed outside the loop, I think

It is pure code and hard to read and to understand. This way I do not see,
why you need Gauss-Hermite.

Could you drop some explanations? For the actual task, if possible, and
the method used

At least both the integrals can be computed outside the loop, I think

Best would be to upload your task as worksheet having the essentials and stating
the task as concrete integral.

Usually there is no need to code Gauss-Hermite, Maple is very good on such things.

Best would be to upload your task as worksheet having the essentials and stating
the task as concrete integral.

Usually there is no need to code Gauss-Hermite, Maple is very good on such things.

Yes, for z < 0 there is no solution, since it would solve N(t)/D(N)(t) * z = sqrt(30),
t = (6*x-z)/sqrt(30), where N = cum. normal and its derivative are positive over IR.

PS: I like it :-)

Yes, for z < 0 there is no solution, since it would solve N(t)/D(N)(t) * z = sqrt(30),
t = (6*x-z)/sqrt(30), where N = cum. normal and its derivative are positive over IR.

PS: I like it :-)

with(Statistics):

dummy := sqrt((1+y*e)/(y*e));
dummy1 := sqrt(1/(1+y*e));
Q := RandomVariable(Normal(0, dummy));
R := RandomVariable(Normal(y*e*z/(1+y*e), dummy1));
y := 1;
e := 1/5; #.2;
k := 1;
#z :=1;
l:= 1/2; #0.5;
a:= 5/4; #1.25;

Int(k*r*PDF(R,r),r=-infinity..x);
#simplify(%);
value(%);
simplify(expand(%), size): combine(%, exp): 
J:=Student[Precalculus]:-CompleteSquare(%,x);

You have to study its zeroes first before proceeding, I would say