Talk:Ziggurat algorithm

Mathematics Start‑class Low‑priority

	Mathematics portal This article is within the scope of WikiProject Mathematics, a collaborative effort to improve the coverage of mathematics on Wikipedia. If you would like to participate, please visit the project page, where you can join the discussion and see a list of open tasks.MathematicsWikipedia:WikiProject MathematicsTemplate:WikiProject Mathematicsmathematics
Start	This article has been rated as Start-class on Wikipedia's content assessment scale.
Low	This article has been rated as Low-priority on the project's priority scale.

This is the talk page for discussing improvements to the Ziggurat algorithm article.
This is not a forum for general discussion of the article's subject.

Put new text under old text. Click here to start a new topic.
New to Wikipedia? Welcome! Learn to edit; get help.

Article policies

Find sources: Google (books · news · scholar · free images · WP refs) · FENS · JSTOR · TWL

Diagram needed

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

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)[reply]

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

The areas A under the curve are equal, and
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:

Choose a point in a vertical interval divided evenly into 8 regions. This gives the slice number i.
Map that region number, via a loojup table, onto a slice of non-uniform height and width.
Choose a point x uniformly between 0 and x_i−1
Test if the point is less than x_i, and accept x immediately if so.
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.
(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)[reply]

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)[reply]