@casperyc You have changed the methodology, and should beware of premature evaluation.

Try this,

restart:

myf := Int(exp(-(x-mu)^2/2/sigma^2)/sigma/sqrt(2*Pi)/(1+exp(-x)),
           x=-infinity..0, epsilon=1e-5, method=_d01amc)
     + Int(exp(-(x-mu)^2/2/sigma^2)/sigma/sqrt(2*Pi)/(1+exp(-x)),
           x=0..infinity, epsilon=1e-5, method=_d01amc):

newmyf:=student[changevar](x=u*sigma,myf,u):

evalf(eval(myf,[mu=-1.645074,sigma=15000])),
evalf(eval(newmyf,[mu=-1.645074,sigma=15000])); # they agree, good

myf0:=unapply('evalf'(myf),[mu,sigma]):
myf0(-1.645074,16500); # does not succeed

newmyf0:=unapply('evalf'(newmyf),[mu,sigma]):
newmyf0(-1.645074,16500); # but this does succeed

plot('newmyf0'(-1.645074,sigma),sigma=17000..18000); # without quotes, you will wait!

plot('newmyf0'(-1.645074,sigma),sigma=1000..18000);

plot3d('newmyf0'(mu,sigma),mu=-10..10,sigma=100..18000,axes=box);

plot3d('newmyf0'(mu,sigma),mu=-10..10,sigma=3000..7000,axes=box);

@casperyc You have changed the methodology, and should beware of premature evaluation.

Try this,

restart:

myf := Int(exp(-(x-mu)^2/2/sigma^2)/sigma/sqrt(2*Pi)/(1+exp(-x)),
           x=-infinity..0, epsilon=1e-5, method=_d01amc)
     + Int(exp(-(x-mu)^2/2/sigma^2)/sigma/sqrt(2*Pi)/(1+exp(-x)),
           x=0..infinity, epsilon=1e-5, method=_d01amc):

newmyf:=student[changevar](x=u*sigma,myf,u):

evalf(eval(myf,[mu=-1.645074,sigma=15000])),
evalf(eval(newmyf,[mu=-1.645074,sigma=15000])); # they agree, good

myf0:=unapply('evalf'(myf),[mu,sigma]):
myf0(-1.645074,16500); # does not succeed

newmyf0:=unapply('evalf'(newmyf),[mu,sigma]):
newmyf0(-1.645074,16500); # but this does succeed

plot('newmyf0'(-1.645074,sigma),sigma=17000..18000); # without quotes, you will wait!

plot('newmyf0'(-1.645074,sigma),sigma=1000..18000);

plot3d('newmyf0'(mu,sigma),mu=-10..10,sigma=100..18000,axes=box);

plot3d('newmyf0'(mu,sigma),mu=-10..10,sigma=3000..7000,axes=box);

@Jarekkk One of the benefits of having Mat_proc_newer accept Matrix Mf and Vector yf as float[8] rtable arguments is that it can be called under evalhf. That's what makes Minimize work fast with it.

And it can also be used that way for single calls, eg,

restart:
NN:=7054:
M := Matrix(NN, 60):
for i to 10 do M[1, i] := 1 end do:
y:=[0.8045, 0.1834, 0.0006]:

Mat_proc_newer := proc(x,M,NN,y)
  local i::integer, j::integer;
  for i from 2 to NN do
    M[i, 1] := M[i-1, 2];
    M[i, 60] := M[i-1, 59];
    for j from 2 to 59 do
      M[i, j] := M[i-1, j] + x * (M[i-1, j-1] - 2*M[i-1, j] + M[i-1, j+1])
    end do
  end do;
return add(((1/10)*add(M[NN, i+10*j], i = 1 .. 10)-y[j])^2, j = 1 .. 3);
end proc:

Mf:=Matrix(M,datatype=float[8]):
yf:=Vector(y,datatype=float[8]):

CodeTools:-Usage(Mat_proc_newer(0.3,Mf,NN,yf));
memory used=79.00MiB, alloc change=11.75MiB, cpu time=1.50s, real time=1.50s
                      0.4419552141658895

CodeTools:-Usage(evalhf(Mat_proc_newer(0.3,Mf,NN,yf)));
memory used=0.58KiB, alloc change=0 bytes, cpu time=265.00ms, real time=263.00ms
                      0.441955214165889509

@Jarekkk One of the benefits of having Mat_proc_newer accept Matrix Mf and Vector yf as float[8] rtable arguments is that it can be called under evalhf. That's what makes Minimize work fast with it.

And it can also be used that way for single calls, eg,

restart:
NN:=7054:
M := Matrix(NN, 60):
for i to 10 do M[1, i] := 1 end do:
y:=[0.8045, 0.1834, 0.0006]:

Mat_proc_newer := proc(x,M,NN,y)
  local i::integer, j::integer;
  for i from 2 to NN do
    M[i, 1] := M[i-1, 2];
    M[i, 60] := M[i-1, 59];
    for j from 2 to 59 do
      M[i, j] := M[i-1, j] + x * (M[i-1, j-1] - 2*M[i-1, j] + M[i-1, j+1])
    end do
  end do;
return add(((1/10)*add(M[NN, i+10*j], i = 1 .. 10)-y[j])^2, j = 1 .. 3);
end proc:

Mf:=Matrix(M,datatype=float[8]):
yf:=Vector(y,datatype=float[8]):

CodeTools:-Usage(Mat_proc_newer(0.3,Mf,NN,yf));
memory used=79.00MiB, alloc change=11.75MiB, cpu time=1.50s, real time=1.50s
                      0.4419552141658895

CodeTools:-Usage(evalhf(Mat_proc_newer(0.3,Mf,NN,yf)));
memory used=0.58KiB, alloc change=0 bytes, cpu time=265.00ms, real time=263.00ms
                      0.441955214165889509

@herclau The original submission has Qnum as an expression, and not a procedure or operator. (I had forgotten that I'd earlier accomodated this using `unapply`, in my answer above.)

So, using Maple 15.01 (64bit, Windows 7), a result better than that from the discrete plot grid is obtained just using,

> Optimization:-Maximize(Qnum, Do = 0.5e-1 .. .5, N = 1 .. 800);

        [24560.2479846227870, [Do = 0.32224998719071896,  N = 142.62810464042423]]

@herclau The original submission has Qnum as an expression, and not a procedure or operator. (I had forgotten that I'd earlier accomodated this using `unapply`, in my answer above.)

So, using Maple 15.01 (64bit, Windows 7), a result better than that from the discrete plot grid is obtained just using,

> Optimization:-Maximize(Qnum, Do = 0.5e-1 .. .5, N = 1 .. 800);

        [24560.2479846227870, [Do = 0.32224998719071896,  N = 142.62810464042423]]

@Thomas Michael Curson We cannot see the embedded image in your post. Is it an image of more code, perhaps pasted (from a Mac?) instead of uploaded using the green up-arrow in the posting editor?

@maple fan It turns out that `int` can also be used here, since it will try numerical intergation if the range values (x below) are floats (or if something like fsolve wraps evalf calls around it, when evaluating at specific points). With Maple 15.01 (Linux, 64bit),

> st:=time():

> fsolve(x->int('g(t)', t=0..x)-5, 0..10);

                                  7.500000000

> time()-st;

                                     0.390

@maple fan It turns out that `int` can also be used here, since it will try numerical intergation if the range values (x below) are floats (or if something like fsolve wraps evalf calls around it, when evaluating at specific points). With Maple 15.01 (Linux, 64bit),

> st:=time():

> fsolve(x->int('g(t)', t=0..x)-5, 0..10);

                                  7.500000000

> time()-st;

                                     0.390

And also see Preben Alsholm's procedure, ResizePlot, here.

And also see Preben Alsholm's procedure, ResizePlot, here.

@hermitian As you've guessed, CodeTools:-Usage is just something that you can wrap around a computational function call, to get an easy printing of the time and memory resources that get used. It's not necessary here, and might not have been available in Maple 13.

The `epsilon` is an option to evalf(Int(...)) the numeric integration (quadrature) command. It is the accuracy tolerance. You see, fsolve is picky and might not return a result unless it is satisified with the convergence of its result. And fsolve might try and attain an accuracy that depends upon the current Digits setting at the level at which it was called. So I set the quadrature accuracy-tolerance `epsilon` to be small enough to get 10 digits of accuracy, but put that inside the EnergyDensity procedure which had Digits raised high enough (hopefully) to allow evalf/Int to satisfy such a tolerance.

The other stuff inside evalf(Int(...)) was to force use of a method that might be suitable for a semi-infinite range of numeric integration. See the evalf/Int help-page. If Digits=15 had not been high enough to allow evalf/Int to get the provided `epsilon` accuracy then I'd have had to set Digits even higher within `EnergyDensity` and not specified _d01amc which only works up to Digits=15 and hardware double precision.

acer 

@hermitian As you've guessed, CodeTools:-Usage is just something that you can wrap around a computational function call, to get an easy printing of the time and memory resources that get used. It's not necessary here, and might not have been available in Maple 13.

The `epsilon` is an option to evalf(Int(...)) the numeric integration (quadrature) command. It is the accuracy tolerance. You see, fsolve is picky and might not return a result unless it is satisified with the convergence of its result. And fsolve might try and attain an accuracy that depends upon the current Digits setting at the level at which it was called. So I set the quadrature accuracy-tolerance `epsilon` to be small enough to get 10 digits of accuracy, but put that inside the EnergyDensity procedure which had Digits raised high enough (hopefully) to allow evalf/Int to satisfy such a tolerance.

The other stuff inside evalf(Int(...)) was to force use of a method that might be suitable for a semi-infinite range of numeric integration. See the evalf/Int help-page. If Digits=15 had not been high enough to allow evalf/Int to get the provided `epsilon` accuracy then I'd have had to set Digits even higher within `EnergyDensity` and not specified _d01amc which only works up to Digits=15 and hardware double precision.

acer 

@tomasz74 If your original statement, "Data points are taken in 1Hz" means that the data points form a grid (regular in both x- and y-directions) then you likely could be able to get a smooth surface from the (stock, non-add-on) plots:-surfdata command. Sorry if that's a repeat, and not appropriate for your situation.

@tomasz74 If your original statement, "Data points are taken in 1Hz" means that the data points form a grid (regular in both x- and y-directions) then you likely could be able to get a smooth surface from the (stock, non-add-on) plots:-surfdata command. Sorry if that's a repeat, and not appropriate for your situation.