@tomleslie 
With a couple of little changes the timing is almost the same:

sinc2:= proc(x) option remember; (sin(x)/x)^2 end:
sinc2(0):=1: sinc2(0.):=1.:
n:=1000;
a:=5.0; h:=a/n;
Z:=Matrix(2*n+1, (i,j) -> sinc2((i-1-n)*h)*sinc2((j-1-n)*h), datatype=float[8],shape=symmetric);

Note that evalhf does not change much because Z has float[8].

The next one is 30% faster:

n:=1000;
a:=5.0; h:=a/n;
sinc3:= proc(i) option remember; local u:=(i-1-n)*h; (sin(u)/u)^2 end:
sinc3(n+1):=1:
Z:=Matrix(2*n+1, (i,j) -> sinc3(i)*sinc3(j), datatype=float[8],shape=symmetric);

@Markiyan Hirnyk 

Cite:

Description
The implicitplot command computes the two-dimensional plot of an implicitly defined curve.  By default, the curve is computed in Cartesian coordinates.

Note that I have added "probably" and "experimental".
That's the "reality".
 

In the help it is mentioned that implicitplot is designed for curves.
It's true that in an example appears a region (f(x,y)<0) but probably this part is only experimental.

@jessica99 

truncation order.

You could have typed ?mtaylor  for Maple help.

Anyway implicitplot should be used only when explicit or parametric are out of the question. So:

plot([seq(arctan(surd(3*cos(x),3))+k*Pi,k=-1..1)], x=-Pi..Pi,y = -(1/2)*Pi-1 .. (1/2)*Pi+1, color=red) ;

@hitstudent 

Maple has it but the user must formulate mathematically the problem.
In this problem Maple will by really helpful if e.g. the direction of the projection is not the normal of the plane or if the sides of the square are not paralel/orthogonal to the base of the cylinder.

@Markiyan Hirnyk 

Thank you. Not a bad job done by Mathematica.
However the two samples N(1,1),  N(3,1/2) seem to be too separated, so the task is easier.

@Carl Love 

I don't think that the command works this way. It could return a convex combination of a normal and some other distribution. It would be interesting to merge two equal samples N(m1,s1), N(m2,s2) and see whether m1,m2,s1,s2 are retrived.

Maybe Markiyan could do this.

@Markiyan Hirnyk 

It was just an opinion, I am not a statistician. However, it does not seem to be  a very complicated command. In your example it returned a convex combination of 2 standard distributions. Having a finite set F of such distributions, it is not difficult to look for a distribution c*f + (1-c)*g,  f,g in F, c in [0,1]  and determine the pair (including parameters) which best fits for a given sample.

@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]);