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*10¹⁷  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  

From the help page

outliers=truefalse  (default)
Outliers are points that are farther than 3/2 times the interquartile range away from the upper and lower quartiles. If this option is set to false, then the outlying points will not be included in the plot. In this case, the whiskers will be extended past outliers to the maximum and minimum points.

In your case
 

restart
with(Statistics):
X := [13, 15, 16, 16, 19, 20, 21, 21, 22, 22, 25, 25, 25, 25, 30, 33, 33, 35, 35, 35, 35, 36, 40, 45, 46, 52, 70]:

IQR := InterquartileRange(X)
                          HFloat(14.0)
LOQ := Quantile(X, 1/4);
UPQ := Quantile(X, 3/4);
                          HFloat(21.0)
                          HFloat(35.0)

OutliersAreIn = [ (-infinity .. LOQ - (3/2)*IQR), (UPQ + (3/2)*IQR .. +infinity) ];
          OutliersAreIn = [-infinity .. HFloat(0.0), HFloat(56.0) .. infinity]

# Thus 70 is considered as an outlier and ths represented by a symbol

The standard boxplot representation is the one used in Maple, see here Box_plot for instance.
As a boxplot is intended to be a graphical summary of a distribution, only its main statistics are represented and its ouliers are treated apart as "rare events".
It is often more informative to see where the outliers are rather than drowning them into the lines.

@Carl Love @Kitonum

Thanks for your explanations.
I should have read the help more carefully.

@acer 

Thank you very much, this will be very useful to me.

@acer 

Thanks you for all these details.

Nevertheless I think I've put you on the wrong track because of my misuse of the term 'indexed" that I should have replaced by "subscripted".
So, I was more or less aware of what you say in your reply.

Now can I rewrite my question this way: 

Is it possible to subscript x by a list?
 

An exemple:
suppose you have the ode x'' + a*x'+b*x=f.
You mightlike to name its solution x__[a, b], and also name x__a the solution which corresponds to b=0

Another example
suppose you have a function f(x, y, z).
You might like to define F__[x, y] as the integral of f over a domain in the [x, y] plane.

I had thouht using MathML to define  the "subscript" (what is right to "__"), but even the simplest 

p := `#mo(a)`;
X__||p

doesn't work
 

@Carl Love 

Hi, 
I often use this syntax and I'm often annoyed to index a symbol not by another symbol but by a couple of symbols.
Here is an example:

• P could represent e pressure,
• T a temperature
• and X some thermodynamic quantity which depends on P and T

This works well

for i in [P, T] do
  X__||i
end do;
                              X__P
                              X__T

but how could I obtain the output

X__[P, T]

Tthe only thing I've been able to do so far uses an intermediate quantity `[P, T]`

h__||`[P, T]`

@vv 

Sorry, I've just read your reply a few seconds after acer's,
I vote up too

@acer 

Thank you acer.
I vote up

@tomleslie 

I think you have missed something:
u is not equal to 

u:= convert(f, string);

but defined by

u  := sprintf("%Zm", f);

and I can't change the way u is constructed.


Where does this sprintf("%Zm", f) come from?
Here is an example

C1 := Cell( Textfield(style=TwoDimOutput,Equation(sqrt(x))) ):
T := Table(Column(), widthmode=percentage, width=10, alignment=center,
           Row( C1)):
InsertContent(Worksheet(Group(Input( T )))):  # Maple 2015

Looking to T, we see that the region corresponding to sqrt(x) corresponds to the code

`_XML_Text-field`("alignment" = "centred", "style" = "2D Output", "layout" = "Normal", _XML_Equation("executable" = "false", "style" = "2D Output", "KiQlInhHIyIiIiIiIw=="))

where 

 "KiQlInhHIyIiIiIiIw=="

is the "code" for sqrt(x).
I was able to find ( showstat(DocumentTools:-Layout:-Equation) ) that this code was obtained this way

u  := sprintf("%Zm", sqrt(x));    

     u := "*$%"xG#"""""#"
               
s  := StringTools:-Encode(u,':-base64');

    s := "KiQlInhHIyIiIiIiIw=="

"Inverting" the last coding to retrieve u is obvious

v := StringTools:-Decode( s, 'encoding' = ':-base64' );

    v := "*$%"xG#"""""#"

Now I wonder if it is possible to "invert" the first coding to retrieve the expression sqrt(x).

Why do I want to do this:

• Suppose someone has constructed the table T and has given you T itself and NOT the code he used to build T.
  in the example above the only thing you know is this

  _XML_Table("interior" = "all", "showinput" = "true", "alignment" = "center", "exterior" = "all", "width" = "10%", "showlabel" = "false", "captionalignment" = "0", "title" = "", "drawtitle" = "false", "order" = "row", "drawcaption" = "false", "randomized" = "false", "captionposition" = "1", "showgroup" = "false", "plotalignlists" = "", "pagebreak" = "none", "postexecute" = "advance", `_XML_Table-Column`("separator" = "true", "weight" = "100"), `_XML_Table-Row`("align" = "top", "separator" = "true", `_XML_Table-Cell`("columnspan" = "1", "backgroundstyle" = "0", "rowspan" = "1", "fillcolor" = "[255,255,255]", `_XML_Text-field`("alignment" = "centred", "style" = "2D Output", "layout" = "Normal", _XML_Equation("executable" = "false", "style" = "2D Output", "KiQlInhHIyIiIiIiIw==")))))

• You can visualise this table and see that the cell contains sqrt(x), but how can you obtain this expression by decoding the las operatot of the field _XML_Equation? 
   
• In this simple cas you could tell me "just write sqrt(x)", but it's not a safe solution for _XML_Equation may contain a more complex expression and the table may contain several _XML_Equation fields.