@TheFixer 

For cones and cylinder "h" is not the same. Look at my vol.

@TheFixer 

In my formula the max volume is 3.175419717

while X=Y=sqrt(2)/2 corresponds to t=Pi/4 and the volume is 2.528158531.

Yes, but probably it would be used only by curiosity, because AFAIK the user usually knows the theoretical distribution of his data and has to determine only the parameters.
Of course it would be interesting to know a statistician's opinion about this.

@nika 

To plot in the same window the histogram and the density function you need relative frequency or scale the density function.

Also, it seems that "a" does not follow very closely a Gamma distribution. (See a true sample in the attached worksheet).

gamma.mw

@nika 

Your data "a" vanished.
It would be a good idea to upload your worksheet.

@J4James 

Actually you don't need, because the fact that f is undefined at those points does not help much:

f:=x->x-floor(x/4)*4:
A:=plot(f, -12..12-0.0001, scaling=constrained, discont):
B:=plot([seq([-12+4*k,0],k=0..5), seq([-8+4*k,4],k=0..5) ],style=point, symbol=circle):
plots[display](A, B, size=[1200,400]);

@torabi 

It seems that you want to obtain some polynomial approximations for an ODE system. You impose e.g. K(0)=0 but in your system K appears at the denominator; it cannot work.

If you change
s2:=K(eta)(0) = 0.001: s3:=Omega(eta)(0) = 0.001:
then fsolve finds a solution, but I don't know if it is related to the solution of the ODE (provided it exists, which is problematic).

@Markiyan Hirnyk 

I don't think that you are in a good position to talk about seriouseness and also politeness here.

@Markiyan Hirnyk 

The post was edited. The Gamma stuff was in the second question of the OP.

@Markiyan Hirnyk 

This is not a scientific paper. The OP assumed that the sample comes from a Gamma distribution and I have shown how the parameters a,b can be obtained. That.s all.

plot/discont seems to be confused when f(x0) is undefined.
This is mathematically acceptable because the continuity makes no sense at a point which is not in the domain of definition.
Actually plot/discont seems to ignore the value of f at x0. Instead, it considers the values near x0, in particular the side limits (if they exist).
For example, f := frac ploted in [0,5] with option discont will have a vertical line at x=5 even if we define f(5):=1 (and so f is now continuous at 5 in [0,5]).

Note that
f:=x->piecewise(x<=4, frac(x), 1);
plots without a vertical line at x=4, ignoring that f(4)=0.

@jamunoz 

Of course the precision depends on the x-step but I don't know about an error estimate in computing the derivative [we would need for example  bounds for f''' ].
We can use higher-order methods for f'(x), e.g.

f'(x)  ~= (f(x-2h) - 4f(x-h) + 3f(x)) / (2h)

@eclark 

It gives more than that, namely the open ends of the interval; in this way tou can replace Open(b) e.g. by b - eps.
You can write a very simple proc to extract a,b and obtain the desired form.

torange:=proc(ineqs::{set,list}, var::name)
  local a,b,x, rr:=subs(var=x,ineqs);
  assume(op(rr));
  rr:=op([1,2],subs[eval]([Open = (a->a)], getassumptions(x)));
  eval(rr, RealRange = ((a,b)-> a..b));
end:

torange({theta <= Pi-arccos(-3/4+(1/4)*sqrt(13+16*sqrt(2))), arccos(3/4) < theta}, theta);

Expand what? The question is more than ambiguous (and your reply too, BTW).

@Matt C Anderson 

@Carl Love 
No, the modules are just fine. As I said, it was pure curiosity, but I don't use them very often (except packages) and I was trying to understand if in such situations they are really more efficient.
Using modules the debugging could be more difficult if something goes wrong.

Edit. You say: "the syntax of the inner procedure is checked when it's originally entered. But the assignment := is not done".
AFAIK, when defined, the body of the procedure is translated into internal code (that is what I understand by parse here). So, why should the inner procedure be "re-parsed"?