Talk:Ziggurat algorithm

Diagram needed

I added an animated diagram.Cuddlyable3 17:58, 14 September 2007 (UTC)

Basically, something like Fig. 1 on p.2 of Marsaglia's paper, or Fig. 5 from the Thomas et al. paper (linked from article References). An improvement would be to add labels "y₀", "y₁", etc. along the y-axis and "x₀", "x₁", etc. above the points.

Without that picture, the best wording is awfully hard to follow.

For a normal distribution, the function to plot is y = e^-x²/2, and the 8-layer Ziggurat points on it are

x[0] = 2.3383716982472527044692   y[0] = 0.0649595117330380202104
x[1] = 1.9819049364005106693290   y[1] = 0.1402998180884519367696
x[2] = 1.7165081257767795707347   y[2] = 0.2291908828832851165061
x[3] = 1.4853586756432938545359   y[3] = 0.3318257830466926345376
x[4] = 1.2629701985308330330774   y[4] = 0.4504325835506035839847
x[5] = 1.0273863717802284625238   y[5] = 0.5899241094185249870813
x[6] = 0.7383689179764494421219   y[6] = 0.7614016031422429793942
x[7] = 0.0000000000000000000000   y[7] = 1.0000000000000000003253

Here's a sample gnuplot input file:

#set terminal svg
#set output "zig.svg"
unset label
set label "x-1" at 2.7120530222393384079012,0.0649595117330380202104 cent offset 0,0.5
set label "x0" at 2.3383716982472527044692,0.1402998180884519367696 cent offset 0,0.5
set label "x1" at 1.9819049364005106693290,0.2291908828832851165061 cent offset 0,0.5
set label "x2" at 1.7165081257767795707347,0.3318257830466926345376 cent offset 0,0.5
set label "x3" at 1.4853586756432938545359,0.4504325835506035839847 cent offset 0,0.5
set label "x4" at 1.2629701985308330330774,0.5899241094185249870813 cent offset 0,0.5
set label "x5" at 1.0273863717802284625238,0.7614016031422429793942 cent offset 0,0.5 
set label "x6" at 0.7383689179764494421219,1.0000000000000000000000 cent offset 0,0.5

set label "y0" at 0,0.0649595117330380202104 right offset -0.5,0
set label "y1" at 0,0.1402998180884519367696 right offset -0.5,0
set label "y2" at 0,0.2291908828832851165061 right offset -0.5,0.2
set label "y3" at 0,0.3318257830466926345376 right offset -0.5,0
set label "y4" at 0,0.4504325835506035839847 right offset -0.5,0
set label "y5" at 0,0.5899241094185249870813 right offset -0.5,-0.7
set label "y6" at 0,0.7614016031422429793942 right offset -0.5,0
set label "y7" at 0,1.0000000000000000000000 right offset -2.5,0

set xrange [0:3.5]
set yrange [0:1.05]
plot exp(-x**2/2) title "exp(-x^2/2)" with lines linewidth 2, "zig.dat" ti "" with lines

And the corresponding input file zig.dat:

0.0000000000000000000000        1.0000000000000000000000
0.7383689179764494421219        1.0000000000000000000000
0.7383689179764494421219        0.5899241094185249870813

0.0000000000000000000000        0.7614016031422429793942
1.0273863717802284625238        0.7614016031422429793942
1.0273863717802284625238        0.4504325835506035839847

0.0000000000000000000000        0.5899241094185249870813
1.2629701985308330330774        0.5899241094185249870813
1.2629701985308330330774        0.3318257830466926345376

0.0000000000000000000000        0.4504325835506035839847
1.4853586756432938545359        0.4504325835506035839847
1.4853586756432938545359        0.2291908828832851165061

0.0000000000000000000000        0.3318257830466926345376
1.7165081257767795707347        0.3318257830466926345376
1.7165081257767795707347        0.1402998180884519367696

0.0000000000000000000000        0.2291908828832851165061
1.9819049364005106693290        0.2291908828832851165061
1.9819049364005106693290        0.0649595117330380202104

0.0000000000000000000000        0.1402998180884519367696
2.3383716982472527044692        0.1402998180884519367696
2.3383716982472527044692        0.0000000000000000000000

0.0000000000000000000000        0.0649595117330380202104
2.7120530222393384079012        0.0649595117330380202104
2.7120530222393384079012        0.0000000000000000000000

71.41.210.146 16:33, 19 June 2007 (UTC)

Thanks for the image, but the caption and animation seem to imply that:

1. The areas A under the curve are equal, and
2. The right-hand (solid white) part is eliminated by rejection.

Neither of those are true. Each layer's black + vertical hatched regions total a constant area A (except for the base layer, which is special), and the right-hand region is eliminated by multiplying a [0,1)-distributed random point by the width of the slice x_i.

I tried to edit the caption to clarify the second point, but the first is pretty hard to fix.

Also, the fact that the distribution tail is, in fact infinite, is not clear from the graphics. It's asymptotic to, but never quite reaches, the X axis.

Sorry to complain, but to illustrate it accurately, you have to demonstrate:

1. Choose a point in a vertical interval divided evenly into 8 regions. This gives the slice number i.
2. Map that region number, via a loojup table, onto a slice of non-uniform height and width.
3. Choose a point x uniformly between 0 and x_i−1
4. Test if the point is less than x_i, and accept x immediately if so.
5. Otherwise, generate a random point y between y_i−1 and y_i and test if y < f(x). If so, accept the point. If not, restart from the beginning.
6. (Step 5 is different in the i=0 case, but let's not try to illustrate that.)

71.41.210.146 02:06, 25 September 2007 (UTC) (originally posted on cuddlyable3's talk page)

I agree that the animation shows

1. The areas A under the curve are equal

That appears to be the case using the 8 coordinates posted here (with high precision!) by anonymous user 71.41.210.146 on 15:26 20 June 2007. The same user now says "Each layer's black and vertical hatched regions total a constant area A" implying x_n(y_n-y_n+1) = constant. That, if true, means his coordinates are unusable. It is impossible to define the base layer where x-->infinity this way, whereas the tail area is finite.

Some other issues in the graphic:

a) The asymptotic tail of the function, as noted. I can remove a white pixel to leave a gap.

b) The left PRNG must select layers with equal probabilities. A non-linear grating suggests what would be done by a look-up table.

c)71.41.210.146 states that each output of the upper PRNG must be multiplied by the x_n of the current layer (so only small wastage) while the graphic shows no conditional multiplication (hence a larger wastage). That seems to trade time to multiply every sample with time to fetch about twice as many PRNs. Cuddlyable3 16:12, 25 September 2007 (UTC)

• It is indeed true that x_n−1(y_n−y_n+1) = constant. The base layer and the tail are handled specially with the i=0 test in step 4, as described in the text leading up to and including the "fallback algorithms for the tail" section. See particularly the description of the "fictitious x₋₁"; the "box" for the base layer is special: shorter than the tail, and of the same area as the base layer under the tail.
• On modern processors, a floating-point multiply is faster than an unpredictable branch, much less a branch plus an additional random number generation (which usually involves a few multiplies itself). Also, when using more than 8 layers, the top layers are very narrow, and rejection would be much more than 50%.

If the special handling of the base layer is not clear from the text, some help improving it would be appreciated! It seems clear enough to me, but if you can point out a place where you get confused, I can try to reword that area.

71.41.210.146 08:41, 29 September 2007 (UTC)

(question removed by Cuddlyable3 17:26, 1 October 2007 (UTC))

Image

The image of of too poor quality. It is better to remove it than keep it in the way it is. Oleg Alexandrov (talk) 01:25, 26 September 2007 (UTC)

Let's talk about how we can make it better. What are your ideas Oleg? Cuddlyable3 10:49, 26 September 2007 (UTC)

• The asymptote of the distribution to the X axis isn't shown well. The distribution extends to infinity.
• The areas marked "A" aren't equal area. You have to include all of the partial box (vertically hatched region) on the non-base layers, specifically including the area above the curve, to get the per-layer area A.
• The way that layer selection is illustrated implies that the probability is proportional to the height. It's not; it's the same for each layer.
• The way the X coordinate is chosen implies that the area to the right of the curve (solid white) is rejected. It's not; the random number is chosen to be uniform from 0 to the width of the layer.

The correct way to draw it in black and white, IMHO, is with solid white boxes under the curve and hollow boxes extending past the curve. Although I don't wish to slight the effort put into the illustration (pictographic illustration of random number generation is extremely challenging!), honestly, a non-animated illustration would do just as well.

But if you like animation, then you have to show that a layer is selected uniformly, and I'd do it like this. Rather than showing a single generation, it shows many random number generations in parallel:

1. Start as shown, with the two-sided distribution cut in half. Excellent, just show the asymptote of the tail.
2. Then show the boxes on top of the curve, and delete the curve everywhere except for the lowest layer.
3. Show that the layers have equal area A.
4. Replace the tail on the lowest layer with a box, preserving the area A.
5. Show a square box on the left, divided horizontally into uniform layers, also each of area A.
6. Fill the square box with random points.
7. Morph the points onto the equal-area layers under the curve. It's just a linear scaling.
8. Extend the right-hand margin of each box down into the box below, dividing it into two pieces (except for the top layer.) Mark the points in each box that are underneath the layer above in green (immediately acceptable). Mark the points in the right-hand portions of each box in yellow (questionable).
9. Make the curve re-appear. On all layers except the bottom, turn yellow points under the curve green. Turn yellow points over the curve red (rejected), then fade them out.
10. Have the yellow points on the bottom layer disappear and reappear as green points in the tail. Don't try to show "morphing" between them.
11. Disappear the boxes and just show the curve and the randomly chosen points under it.

That actually depicts the algorithm quite faithfully.

71.41.210.146 08:26, 29 September 2007 (UTC)

Thank you for your script for a new animation. The .gif animation format is poorly equipped for morphing, and with a bigger palette for colours and fading the file size will swell.Cuddlyable3 16:15, 1 October 2007 (UTC)