You seem to be building an elaborate structure for testing the reproducibility
of the random number generator. I suspect that rbinom is calling the random
number generator a different number of times when you pass prob=0.5 than
otherwise.
---------------------------------------------------------------------------
Jeff Newmiller The ..... ..... Go Live...
DCN:<jdnewmil at dcn.davis.ca.us> Basics: ##.#. ##.#. Live
Go...
Live: OO#.. Dead: OO#.. Playing
Research Engineer (Solar/Batteries O.O#. #.O#. with
/Software/Embedded Controllers) .OO#. .OO#. rocks...1k
---------------------------------------------------------------------------
Sent from my phone. Please excuse my brevity.
Michael Hannon <jm_hannon at yahoo.com> wrote:
>Greetings.? My wife is teaching an introductory stat class at UC
>Davis.? The
>class emphasizes the use of simulations, rather than mathematics, to
>get
>insight into statistics, and R is the mandated tool.?? A student in the
>class
>recently inquired about different approaches to sampling from a
>binomial
>distribution.? I've appended some code that exhibits the idea, the gist
>of
>which is that using sample(c(0, 1), ...) and rbinom(...) should give
>equivalent results.
>
>The surprising (to me) result is that the two approaches DO give the
>same
>result, EXCEPT when the probability is exactly 0.5.? See Appendix A for
>the
>code and Appendix B for the output.? I don't think this issue is
>system-dependent, but I've put my session information in Appendix C.
>
>Another wrinkle in this is that if I omit the "prob" parameter
from the
>call
>to sample, meaning to take the default value of 0.5, the two methods DO
>give
>the same result.
>
>Any thoughts about this?? Thanks.
>
>--Mike
>
>Appendix A: some R code that exhibits the problem
>================================================>
>ppp <- seq(0, 1, by = 0.01)
>
>result <- do.call(rbind, lapply(ppp, function(p) {
>??? set.seed(1)
>??? sampleRes <- sample(c(0, 1), size = 1, replace = TRUE,
>??????????????????????? prob=c(1-p, p))
>?? ?
>??? set.seed(1)
>??? rbinomRes <- rbinom(1, size = 1, prob = p)
>?? ?
>??? data.frame(prob = p, equivalent = all(sampleRes == rbinomRes))
>?? ?
>}))
>
>result
>
>
>Appendix B: the output from the R code
>=====================================>
>??? prob equivalent
>1?? 0.00?????? TRUE
>2?? 0.01?????? TRUE
>3?? 0.02?????? TRUE
>4?? 0.03?????? TRUE
>5?? 0.04?????? TRUE
>6?? 0.05?????? TRUE
>7?? 0.06?????? TRUE
>8?? 0.07?????? TRUE
>9?? 0.08?????? TRUE
>10? 0.09?????? TRUE
>11? 0.10?????? TRUE
>12? 0.11?????? TRUE
>13? 0.12?????? TRUE
>14? 0.13?????? TRUE
>15? 0.14?????? TRUE
>16? 0.15?????? TRUE
>17? 0.16?????? TRUE
>18? 0.17?????? TRUE
>19? 0.18?????? TRUE
>20? 0.19?????? TRUE
>21? 0.20?????? TRUE
>22? 0.21?????? TRUE
>23? 0.22?????? TRUE
>24? 0.23?????? TRUE
>25? 0.24?????? TRUE
>26? 0.25?????? TRUE
>27? 0.26?????? TRUE
>28? 0.27?????? TRUE
>29? 0.28?????? TRUE
>30? 0.29?????? TRUE
>31? 0.30?????? TRUE
>32? 0.31?????? TRUE
>33? 0.32?????? TRUE
>34? 0.33?????? TRUE
>35? 0.34?????? TRUE
>36? 0.35?????? TRUE
>37? 0.36?????? TRUE
>38? 0.37?????? TRUE
>39? 0.38?????? TRUE
>40? 0.39?????? TRUE
>41? 0.40?????? TRUE
>42? 0.41?????? TRUE
>43? 0.42?????? TRUE
>44? 0.43?????? TRUE
>45? 0.44?????? TRUE
>46? 0.45?????? TRUE
>47? 0.46?????? TRUE
>48? 0.47?????? TRUE
>49? 0.48?????? TRUE
>50? 0.49?????? TRUE
>51? 0.50????? FALSE
>52? 0.51?????? TRUE
>53? 0.52?????? TRUE
>54? 0.53?????? TRUE
>55? 0.54?????? TRUE
>56? 0.55?????? TRUE
>57? 0.56?????? TRUE
>58? 0.57?????? TRUE
>59? 0.58?????? TRUE
>60? 0.59?????? TRUE
>61? 0.60?????? TRUE
>62? 0.61?????? TRUE
>63? 0.62?????? TRUE
>64? 0.63?????? TRUE
>65? 0.64?????? TRUE
>66? 0.65?????? TRUE
>67? 0.66?????? TRUE
>68? 0.67?????? TRUE
>69? 0.68?????? TRUE
>70? 0.69?????? TRUE
>71? 0.70?????? TRUE
>72? 0.71?????? TRUE
>73? 0.72?????? TRUE
>74? 0.73?????? TRUE
>75? 0.74?????? TRUE
>76? 0.75?????? TRUE
>77? 0.76?????? TRUE
>78? 0.77?????? TRUE
>79? 0.78?????? TRUE
>80? 0.79?????? TRUE
>81? 0.80?????? TRUE
>82? 0.81?????? TRUE
>83? 0.82?????? TRUE
>84? 0.83?????? TRUE
>85? 0.84?????? TRUE
>86? 0.85?????? TRUE
>87? 0.86?????? TRUE
>88? 0.87?????? TRUE
>89? 0.88?????? TRUE
>90? 0.89?????? TRUE
>91? 0.90?????? TRUE
>92? 0.91?????? TRUE
>93? 0.92?????? TRUE
>94? 0.93?????? TRUE
>95? 0.94?????? TRUE
>96? 0.95?????? TRUE
>97? 0.96?????? TRUE
>98? 0.97?????? TRUE
>99? 0.98?????? TRUE
>100 0.99?????? TRUE
>101 1.00?????? TRUE
>
>Appendix C: Session information
>==============================>
>> sessionInfo()
>R version 3.0.0 (2013-04-03)
>Platform: x86_64-redhat-linux-gnu (64-bit)
>
>locale:
>?[1] LC_CTYPE=en_US.UTF-8?????? LC_NUMERIC=C???????????? ?
>?[3] LC_TIME=en_US.UTF-8??????? LC_COLLATE=en_US.UTF-8?? ?
>?[5] LC_MONETARY=en_US.UTF-8??? LC_MESSAGES=en_US.UTF-8? ?
>?[7] LC_PAPER=C???????????????? LC_NAME=C??????????????? ?
>?[9] LC_ADDRESS=C?????????????? LC_TELEPHONE=C?????????? ?
>[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C????? ?
>
>attached base packages:
>[1] stats???? graphics? grDevices utils???? datasets? methods?? base???
>?
>>
>
>______________________________________________
>R-help at r-project.org mailing list
>https://stat.ethz.ch/mailman/listinfo/r-help
>PLEASE do read the posting guide
>http://www.R-project.org/posting-guide.html
>and provide commented, minimal, self-contained, reproducible code.