{"id":753,"date":"2017-11-26T10:32:45","date_gmt":"2017-11-26T15:32:45","guid":{"rendered":"http:\/\/www.clayford.net\/statistics\/?p=753"},"modified":"2020-09-13T10:17:55","modified_gmt":"2020-09-13T14:17:55","slug":"recreating-a-geometric-cdf-plot-from-casella-and-berger","status":"publish","type":"post","link":"https:\/\/www.clayford.net\/statistics\/recreating-a-geometric-cdf-plot-from-casella-and-berger\/","title":{"rendered":"Recreating a Geometric CDF plot from Casella and Berger"},"content":{"rendered":"

For reasons I no longer remember I decided a while back to reproduce Figure 1.5.2 from Casella and Berger<\/a> (page 32):<\/p>\n

\"\"\/<\/p>\n

It's a plot of the cumulative distribution function of a geometric distribution with p = 0.3<\/em>. The geometric distribution is best explained with coin-flipping. Let's say we keep flipping a coin until we observe a head. Now let X<\/em> be a random variable that equals the number of tosses required to see a head. So if we get a head on the first flip of the coin, that means we needed one toss. If it takes two tosses, we got a tail on the first toss and then a head on the second toss. The graph above visualizes the possibilities for 15 flips with probability of getting heads set to 0.3. We see that the probability of taking only one toss is 0.3, which makes sense since the probability of getting heads is 0.3. The probability of requiring two tosses is about 0.5. And so on. The straight lines indicate the probabilities only change at the whole numbers. There is no such thing as, say, 2.23 flips. This is often referred to as a step function. We see that there's a high probability of getting a head by the 5th or 6th flip. This is a simple waiting-time distribution that can be used to model the number of successes before a failure and vice-versa.<\/p>\n

To recreate this plot in R I used the pgeom<\/em> function, which returns the cumulative probability of the geometric distribution for a given number of failures before a success occurs and a given probability of success. For example, the probability of 0 failures before observing a success with p = 0.3<\/em>:<\/p>\n

pgeom(0, prob = 0.3)\r\n<\/code><\/pre>\n
## [1] 0.3\r\n<\/pre>\n
The probability of 1 failure before observing a success with p = 0.3<\/em>:<\/p>\n
pgeom(1, prob = 0.3)\r\n<\/code><\/pre>\n
## [1] 0.51\r\n<\/pre>\n
To generate all the probabilities in the plot, we can do the following:<\/p>\n
y <- pgeom(q = 0:14, prob = 0.3)\r\n<\/code><\/pre>\n
Notice that instead of using 1:15<\/code> (number of total flips for success), we use 0:14<\/code>. That's because pgeom<\/em> works with the number of failures before<\/em> success. Now that we have our y coordinates, we can create our first version of the plot as follows:<\/p>\n
plot(x = 1:15, y, pch = 19)\r\n<\/code><\/pre>\n
\"plot<\/p>\n
Now let's add the horizontal lines using the segments<\/code> function:<\/p>\n
plot(x = 1:15, y, pch = 19)\r\nsegments(x0 = 0:15, y0 = c(0, y),\r\n         x1 = 1:16, y1 = c(0, y), lwd = 2)\r\n<\/code><\/pre>\n
\"plot<\/p>\n
The x0<\/code> and y0<\/code> coordinates are where the line starts. The x1<\/code> and y1<\/code> coordinates are where the line ends. Since the lines are horizontal, the y coordinates are the same for the start and end postions. The y coordinates include 0, so we add that value to y<\/code> with the c<\/code> function.  The lwd = 2<\/code> argument makes the line a little thicker and darker.<\/p>\n
Our plot is basically done, but just for fun I wanted to see how close I could make it look like the version in the book. That means relabeling the axes, moving the axis labels to the ends, and removing the lines at the top and right side of the plot. It also means moving the axis tick marks inside<\/em> the plotting area. After some trial and error and deep dives into R documentation, here's what I was able to come up with:<\/p>\n
plot(x = 1:15, y, pch = 19, \r\n     yaxt = "n", xaxt = "n", ylim = c(0,1.05), xlim = c(0,15.5), \r\n     bty="l", xaxs="i", yaxs = "i", xlab = "", ylab = "")\r\nsegments(x0 = 0:15, y0 = c(0, y),\r\n         x1 = 1:16, y1 = c(0, y), lwd = 2)\r\naxis(side = 1, at = 0:15, \r\n     labels = 0:15, tcl = 0.5, family = "serif")\r\naxis(side = 2, at = seq(0,1,0.1), \r\n     labels = c(0,paste(".",1:9),1), las=1, tcl = 0.5, family = "serif")\r\nmtext(text = expression(italic(x)), side = 4, \r\n      at = 0, las = 1, line = 0.5, family = "serif")\r\nmtext(text = expression(italic(F[x](x))), side = 3, \r\n      at = 0, line = 0.5, family = "serif")\r\n<\/code><\/pre>\n
\"plot<\/p>\n
In the plot<\/code> function: <\/p>\n