acer

32348 Reputation

29 Badges

19 years, 329 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@qilianshan The image attached above by the Original Poster almost makes it look like he's only interested in the grid of points in the domain consisting of the end-points and mid-points. Obviously there is no need to use a point-probe to get those coordinates, since they can be immediately obtained from the plot command's own input ranges.

So it seems reasonable to conclude that the OP is asking about point-probe of a 3D surface plot.

Somewhere on a harddrive I have a procedure that accepts a plot3d structure and embeds an app with tools for doing that. But before I can post it I have to find it, and that won't be for at least a day.

[edit] FYI, what I remember of the applet I wrote is that it consists of a procedure which you call with the 3D plot as argument. Doing so would automagically display the 3D plot as usual, with a small button to reveal or hide its point-probe tools. That consisted of a 2D plot with the right click&drag trackers, and an equivalent densityplot background from the 3D plot data. And subsequent sliding around of the point-probe would reveal the x-y point's interpolated data as well as overlay a moving marker point and line/transparent-plane on the 3D and 2D plots. But it's on the harddrive of a broken computer in my basement, and it's a beautiful sunny day for a hike with my kids.

@Christopher2222 Perhaps your main irritation is that you don't want to have to move the cursor to the start of the name, and then back again, if you're in the middle of typing in a name and suddenly come to a point where you suddenly realize that you need some unusual characters. In that situation you can use || to concatenate. Eg,

ini||`.`||param;
                        ini.param

lprint(%);
 `ini.param`

Some other ways, that don't involve typing in the left-quote (though not particularly terse):

ini||"."||param;
                        ini.param

lprint(%);
 `ini.param`

convert("ini.param", name);
                        ini.param

lprint(%);
 `ini.param`

nprintf("%a.%a",ini,param);
                        ini.param

lprint(%);
 `ini.param`

nprintf("ini.param");
                        ini.param

lprint(%);
 `ini.param`

@Christopher2222 It really ought to be obvious to you that some mechanism is necessary to disambiguate between the name you want and the multiplication.

@tomleslie Thanks for shedding valuable light on the situation, Tom.

@Carl Love The hypergeom appears to be fast and robust, with the `p` value converted to exact rational before passing to NextZero.

Here is a procedure for it, and some timings, and a loglog plot as Carl showed.
 

restart;

 

Quant:=proc(r,p,P)
  local ans,func,i,pp,res,special,st,t;
  uses RF=RootFinding;
  special:=1-GAMMA(r+t+1)*pp^r*(1-pp)^(t+1)
             *hypergeom([1,r+t+1],[t+2],1-pp)/(GAMMA(r)*GAMMA(t+2));
  Digits:=max(16,Digits);
  UseHardwareFloats:=false;
  # Find value to start NextZero's search
  for i from 1 to 40 do
    res:=evalf(eval(special,[pp=p,t=2^i])-P);
    if res>0 then
      st:=i; break;
    end if;
  end do:
  func:=unapply(simplify(eval(special,
                              [pp=convert(p,rational,exact)])),t)
      assuming t>1:
  ans:=ceil(RF:-NextZero(t->func(t)-P,
                         2^(st-1),
                         ':-maxdistance'=2^st-2^(st-1)));
end proc:

 

CodeTools:-Usage( Quant(2, 3e-3, 0.98) );

memory used=22.22MiB, alloc change=34.00MiB, cpu time=186.00ms, real time=186.00ms, gc time=7.58ms

1941

CodeTools:-Usage( Quant(1.965, 3e-3, 0.98) );

memory used=45.02MiB, alloc change=4.00MiB, cpu time=319.00ms, real time=319.00ms, gc time=22.00ms

1920

CodeTools:-Usage( Quant(1.965, 3e-10, 0.98) );

memory used=54.88MiB, alloc change=-2.00MiB, cpu time=426.00ms, real time=426.00ms, gc time=56.71ms

19239732630

 

plots:-loglogplot(p->Quant(1.965, p, 0.98), 1e-7..0.5,
                  labels = [p, `98th %tile`],
                  labeldirections = [horizontal,vertical],
                  title = "The 98th %tile of NegativeBinomial(1.965, p) vs p",
                  adaptive=false, numpoints=49,
                  smartview=false, gridlines=false);

Download Quant01.mw

@mmcdara I don't understand your point. You used 2 instead of 1.965.

Using 2 instead of 1.965 then the regular Statistics package find the same results quickly. And using 1.965 instead of 2 and the Student:-Statistics package also take a long time...

restart:
with(Statistics):
X := RandomVariable(NegativeBinomial(2, 0.003)):

Quantile(X,0.98, numeric);
                             1941.

X := RandomVariable(NegativeBinomial(2, 3e-10)):
Quantile(X,0.98, numeric);
                                        10
                     1.94464056686655 10  

restart:
with(Student:-Statistics):
X := NegativeBinomialRandomVariable(1.965, 0.003):

Quantile(X,0.98, numeric); # goes away...

 

@David1 The GlobalSearch command searches for local and global extrema, and hopefully finds global extrema. If you pass it the option solutions=1 then it may exit as soon as it finds one. I can imagine that that would decrease the chance of its finding a global extremum. A reasonable interpretation of the help page for this command is that passing solutions=1 could allow it to return as soon as it found a single local extreme point, which is certainly not at all your goal.

If you want DirectSearch to try and find a single point which is the global extremum then why not use its command GlobalOptima?  That command's goal is to search for the best objective value. As you will see on its help page that command takes many options, because (as has been pointed out to you before) there is no foolproof method for computing the globally best result of a general, nonconvex, nonlinear (and possibly discontinuous) optimization problem.

If you want to increase the chances that it finds the global extremum then make it work harder, searching from more randomized starting points and with a longer time-limit and evaluation-limit. That relates to its options number, timelimit, and evaluationlimit. But other options can also come into play. Read the help page carefully. Experiment with the options.

@jamalator How can y and q both be prime and satisfy 8695*y = 1341*q ?

@Joe Riel Since this discussion relates to thismodule then perhaps it's on-topic to mention function:-somestaticexport. (or the subtle issues in invoking :-eval(self) or mapping a commonly named static across different objects in a list -- including doing so in tricky situations).

Sorry, Carl. That was me. Unfortunately there's no way at present for me to toggle the Product entry to a specific Maple version without Mapleprimes considering the Question as being edited. I figured that if I did them in a block it might be more clear that it was a series of no-ops. I try not to do that sort of edit often.

@Melvin Brown You're welcome, Marvin.

@waseem What do you want to combine? With respect to what do you want to collect? Why are you not being clear?

Perhaps you are looking for something like this form?

restart:

with(DETools):
H:=proc(ee)
     local temp;
     temp:=sort([op(indets(ee,And(polynom(integer,r),
                                  satisfies(u->degree(u,r)>0))))],
                (a,b)->degree(a,r)>degree(b,r) and length(a)>length(b));
     subs(map(u->u=freeze(u),temp),ee);
end proc:

P[o](z):=C[o]*exp(lambda*z)+D[o]*exp(-lambda*z):
u[o](r,z):=(1-r^2)*diff(P[o](z),z):
v[o](r,z):=(2*r-r^3)*diff(P[o](z),z,z):

A1:=P[o](z)*u[o](r,z):
A2:=lambda^2*P[o](z)*(1+r^2)+v[o](r,z):

v[1](r,z):=(gamma1^2/16)*diff(P[o](z),z,z)*(2*r-r^3):

u[1](r,z):=G*lambda*(C1^2*exp(2*lambda*z)-D1^2*exp(-2*lambda*z))*(((1-r^4)/16)-((1-r^6)/72))+(lambda/2)*(C1*exp(lambda*z)-D1*exp(-lambda*z))*(1-r^4)-A1*(1-r^2)/4:

gamma1:=4/lambda:
P1:=gamma1^2*diff(u[o](r,z),z,z)+diff(u[1](r,z),r,r)+(1/r)*diff(u[1](r,z),r)-G*(gamma1^2*u[o](r,z)*diff(u[o](r,z),z)+v[1](r,z)*diff(u[o](r,z),r)):

P2:=int(P1,z);

((1/12)*G*lambda*(r^5-3*r^3)*((1/2)*C1^2*exp(2*lambda*z)/lambda+(1/2)*D1^2*exp(-2*lambda*z)/lambda)+(-r^2+1)*r*((1/2)*(exp(lambda*z))^2*C[o]^2+(1/2)*D[o]^2/(exp(lambda*z))^2)-2*lambda*r^3*(C1*exp(lambda*z)/lambda+D1*exp(-lambda*z)/lambda))/r+(-r^2+1)*((1/2)*(exp(lambda*z))^2*C[o]^2+(1/2)*D[o]^2/(exp(lambda*z))^2)+(1/12)*G*lambda*(5*r^4-9*r^2)*((1/2)*C1^2*exp(2*lambda*z)/lambda+(1/2)*D1^2*exp(-2*lambda*z)/lambda)-2*G*(8*(-r^2+1)^2*((1/2)*C[o]^2*lambda^2*(exp(lambda*z))^2+(1/2)*D[o]^2*lambda^2/(exp(lambda*z))^2)/lambda^2-(-r^3+2*r)*r*((1/2)*C[o]^2*lambda^2*(exp(lambda*z))^2+(1/2)*D[o]^2*lambda^2/(exp(lambda*z))^2)/lambda^2)+16*(-r^2+1)*(C[o]*lambda^2*exp(lambda*z)+D[o]*lambda^2*exp(-lambda*z))/lambda^2-6*lambda*r^2*(C1*exp(lambda*z)/lambda+D1*exp(-lambda*z)/lambda)-2*r^2*((1/2)*(exp(lambda*z))^2*C[o]^2+(1/2)*D[o]^2/(exp(lambda*z))^2)

#P[1](z):=simplify(P2);

#ans2:=thaw(factor(combine(expand(H(P2)))));

temp3:=indets(P2,And(polynom(anything,{C1,C[o]}),`+`,satisfies(u->not has(u,r)))):
Rule3:=map(u->u=(uu->K(freeze(numer(uu)))/denom(uu))(factor(combine(expand(u)))),temp3):
ans3:=thaw(eval(collect(subs(Rule3,P2),K,simplify),[K=(()->args)]));

(-2*r^2+1)*(exp(2*lambda*z)*C[o]^2+D[o]^2*exp(-2*lambda*z))-8*r^2*(C1*exp(lambda*z)+D1*exp(-lambda*z))+(1/4)*G*r^2*(r^2-2)*(C1^2*exp(2*lambda*z)+D1^2*exp(-2*lambda*z))+(-16*r^2+16)*(C[o]*exp(lambda*z)+D[o]*exp(-lambda*z))-G*(9*r^4-18*r^2+8)*(exp(2*lambda*z)*C[o]^2+D[o]^2*exp(-2*lambda*z))

simplify(P2-ans3);

0

 

Download waseem_2.mw

If this was your actual goal then why didn't you say so in your original question?

@Carl Love I didn't have a particularly strong reason, except as a single sanity check. (I'm sometimes a little hesitant for trig integrals over a symmetric range about zero, slightly more so when I see the multiple of Pi on the argument.)

I see that Rouben found the same thing, while I was working on it and submitting.

Yet another (much slower) approach that works here is to change variables to make the inner integral range from a=1..infinity, and then do purely numeric integration. Mariusz's nice workaround also serves as another corroboration.

restart;
with(IntegrationTools):
H:=Int(exp(-a^2-b^2-c^2),
       [a=b^2/(4*c)..infinity, b=-infinity..infinity, c=0..infinity]):
T:=ExpandMultiple(H):
T1:=Change(op([1,1],T),a=s*(b^2/(4*c)),[s])
    assuming c>0, b::real:
new:=CollapseNested(Int(Int(T1,op([1,2],T)),op(2..,T))):
evalf(op(0,new)(op(new),epsilon=1e-5));
                          0.9786008269
First 240 241 242 243 244 245 246 Last Page 242 of 592