mmcdara

1781 Reputation

14 Badges

4 years, 233 days

MaplePrimes Activity


These are replies submitted by mmcdara

@acer 

So it was a well known crash cause.
I'm happy to read that it is corrected (even if, in the example I posted, I believe that  is better to avoid doing this stupid thing of sorting a null vector)

Are these crashes recorded somewhere?

They are really borring when you meet one, specially if they occur during the execution of a Maplet-based application (debuggibg is impossible AFAK)

Thanks for your comment

@vv 

"... the maximum point found was outside the interval of integration"
You're right I'd missed this point.
Thanks for pointing it out



 

@acer 

Thank you acer. 
Being iused to end all my lines with a colon or a semicolon (even if it's not necessary), I regularly get trapped by putting an after proc(...)

@vv 

Great, I hadn't had the idea to represent f by a piecewise function, nor had I even thought that it would have give a general solution.

Thanks

PS: if you have, let's say 10 thousand points and "only" 10 "classes" (here the second element od any List[i]), it's better to display 10 plots where each correspond to a points in the same "class", than to display 10 thousands points.

Here is an example with 12 classes
 

restart:

with(ColorTools):
with(plots):
with(Statistics):

r := rand(0.0 .. 1.0):
N := 10000:
X := Vector(N, i->r()):
Y := Vector(N, i->r()):
Z := X +~ Y: # +~ Vector(N, i->r()):

# these 2 lines mimic what a simple colorscheme does
c := EvenSpread(Color("RGB", "Yellow"), 12):
c := op~(1, NearestNamedColor~(c0, palette="X11"));

C   := floor~(Z *~ 12) +~ 1:
eps := (max-min)(Z) * 1e-8:
mZ  := min(Z)-eps:
MZ  := max(Z)+eps:
C := floor~( (Z -~ mZ) /~ (MZ-mZ) *~ 12) +~ 1:


# this is extremely lengthy

CodeTools:-Usage(
  display( seq(plot([[X[n], Y[n]]], style=point, symbol=solidcircle, symbolsize=15, color=c[C[n]]), n=1..N) )
);

LC := Vector(12, i->NULL):
for n from 1 to N do
  LC[C[n]] := LC[C[n]] , n
end do:

# while this is very fast
CodeTools:-Usage(
       display(
         PLOT(
           seq(
             POINTS(
                < X[[LC[n]]] | Y[[LC[n]]]>, 
                COLOR(RGB, op(NameToRGB24(c[n])/~255))
             ),
             n=1..12
           )
         )
       )
     ):

);

memory used=10.19MiB, alloc change=0 bytes, cpu time=80.00ms, real time=117.00ms, gc time=0ns

 

@acer 

Thanks.

Regarding your last remark, I noticed that a lot of questions invoking Maplets remain unanswered. Probably because some think of them as an old-fashioned or exotic technology?
You may not think this is quite honest, but I thought that omitting the word "Maplets" would increase the probability of getting an answer.

 

@vv 

Thank you vv.
I'm completely ignorant of all this, so I'm going to look to to the Autolt Script Editor to see it can be of any help

@acer 

OK, thank you acer.

The specific topic I'm working on is the development of a (simplified) post processor to help people to analyze statistical data.
A strong constraint is that these people have only basic knowledge in Maple.

To do this I used the Maplets technology (which seems more intuitive to the users than Embeded Components).
At one point a 3D plot is displayed in a  Maplets[Plotter] element. The user can rotate it up to an informative point of view and then decide (by clicking on a specific Maplet[Button] element) to export it.

The problem is that this is not the displayed graph which is exported but the one I originally displayed in the  Maplets[Plotter] element. More of this it doesn't seem possible to get the options of the current displayed graph (I can only get what I sent to Maplets[Plotter] ).


It wouldn't matter if I had displayed the graph in the spreadsheet because a right click on the mouse would then allow me to export the displayed graph. But no access to the options of the plot are avaliable in a  Maplets[Plotter] region.
The solution I use for the moment is to add three  Maplets[Slider] elements to make varying theta, phi and psi: this allows the iuser to rotate the figure and their values can be captured (with Maplets[Tools][Get]) and used to build an orientation[...] option I add to the graph before exporting it.
Unfortunately rotating a graph by moving these three Maplets[Slider] is far less natural than to rotate it manually.


About your last phrase "I have requested that Property for Plot Components, some years ago. It might help if someone else who wanted it were to also submit the request." I agree one hundred percent.
Could you tell me how I can do wuch a request? Through a simple question or by an other mean?

 

I found this for you 
to be read Periodization_discretization.html  and   matdid577306.pdf
read it if you want Dirac_comb

@Umang Varshney 

I hadn't looked closely at your definition of d: it contains a BIG mistake.
You write

restart:
with(Statistics):
d := Distribution(
       PDF = (proc (t) options operator, arrow; f__D(t) end proc), 
       Quantile = (q -> int(f__D(t), t = -infinity .. q))
    ):

x := RandomVariable(d);

But a Quantile of order q is not the integral of the PDF from -oo to q:

# For continuous pdf, z is the Quantile of order q of x 
# if z is the only value such that 
int(f__D(t), t = -infinity .. z)=q:

# thus:
Quantile(x, q) = solve(int(f__D(t), t = -infinity .. z)=q, z)

# and a correct definition of d is
#( the ''...'' are necessary to prevent MAPLE to try solving)

d := Distribution(
       PDF = (proc (t) options operator, arrow; f__D(t) end proc), 
       Quantile = (q -> ''solve(Int(f__D(t), t = -infinity .. z)=q, z)'')
     ):
x := RandomVariable(d);

# Quantile of order p (p is a probability)

Quantile(x, p);
                  /  /z                        \
                  | |                          |
             solve| |         f__D(t) dt = p, z|
                  | |                          |
                  \/-infinity                  /

What you have defined with 

Quantile = (q -> int(f__D(t), t = -infinity .. q))

is the Probability that x be less or equal to q.

I think nothing is very serious because I have the impression that you simply used the term "Quantile" instead of the term "Probability".


So my advice is to redefine d (and beta) this way

restart:
with(Statistics):
d := Distribution(
       PDF = (proc (t) options operator, arrow; f__D(t) end proc), 
       CDF = (q -> int(f__D(t), t = -infinity .. q))
    ):

x := RandomVariable(d);


# and latter in my code ABCDE.mw:

F := t -> eval(CDF(x, s), s=t):
F__c := t -> 1- F(t);   # which seems consistent for 0 <= F(t) <= 1


# idem deeper for G(t)


Here is my code corrected

ABCDEF.mw

@Carl Love 

I agree, Quantile (in d and d1) are not correctly defined and the piecewise definition of Pi__ER can be simplified.
Concerning the computation time some rewritting seems necessary.
But, at the very end, as  @Axel Vogt quoted, the question looks like more to the description of a general methodology to be applied to a certain class of problems than to an actual problem, e.g.:

  1. calculate the mean of Pi__ER
  2. find q__ER 
    TIP: derive the mean according to q  and find q__ER such that derivative is null
  3. plug q__ER into the expression of the mean of Pi__ER

 

@Carl Love 

You write "This procedure does the same as my Answer above": yes indeed, I reused your key idea (IMO) about the reduction of n by a factor 2^(k-1) to spare time.
Next Maple 2015 does not have a procedure equivalent to SumOfDivisors, so the remaining part of my code

combinat:-powerset(F):
add(map(u -> mul(op(u)), %))*k

To conclude I do not pretend to have done something original.

Concerning your second reply "Your timing is not accurate because CodeTools:-Usage doesn't account for the fact that the result of ifactor (and thus also ifactors, essentially) is cached": this is very interesting. I have already be surprised by the 
results  CodeTools:-Usage and I generally use time() instead before and after a loop were the code to test is executed many times.
Maybe this is also due to this notion of "cache".
Is this mentionned somewhere in the help pages?
Interesting remark 

@nm @Axel Vogt  

None of your solutions is correct:

  • nm: your solution cannot be plotted for a larger range than the one you used
    (probably numerical problems when evaluating it)
  • Axel: your solution doesn't verify the asymptotic regime (outside of the boundary layer -0.065..0.05)

Here is an approximate solution based on spline approximation of f in the "boundary layer" -0.065..0.05

integral.mw

@Carl Love 

Understood, sorry.


By the way, as I can't use your coe with my "old" Maple 2015 version, I've written this one.
Almost as fast as yours , but the memory size is significanly larger for mine.
 

restart
interface(version)

Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895


F := proc(n)
 local f, k, F:
 f := ifactors(n)[2];
 if f[1][1]=2 then 
   k := `^`(op(f[1])): 
   F := map(u -> `$`(op(u)), f[2...-1])
 else 
   k := 1: 
   F := map(u -> `$`(op(u)), f)
 end if:
 combinat:-powerset(F):
 add(map(u -> mul(op(u)), %))*k
end proc:

n := 2^39*3*5*7*11*13*17*19*23*nextprime(2^29):

CodeTools:-Usage( F(n), iterations=1000 );

memory used=169.74KiB, alloc change=0 bytes, cpu time=967.00us, real time=970.00us, gc time=94.25us

                 82255314605128880924739502080

 

@Carl Love 

Hi Carl,
I'm surprised by your result

restart
interface(version)
Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895

with(numtheory):

n:= 2^39*3*5*7*11*13*17*19*23*nextprime(2^29)

                 32922697295031156088441405440

add(select(is, divisors(n), odd));
                       149621545652060160

n := 3*5*7*11*13*17*19*23*nextprime(2^29);
add(divisors(n));  

                       59886037515809505
                       149621545652060160


 

The sum of all the divisors is about 3.9*1017  thus one can bet the sum of all odd divisors is about half this value

add(divisors(n)):
evalf(%);
                      
                                      17
                        1.496215457 10  

# Gronwall (1913) proved this bound 
# limsup( add(divisors(n)) / n*log(log(n))) = exp(0.5772156649)

evalf( n*log(log(n))*exp(0.5772156649) )
                                      17
                        3.897471248 10  

 

 

 

1 2 3 4 5 6 7 Last Page 1 of 62