<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/mysteriously-slow-sample/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://meanmean.me/2018/07/28/Mysteriously-slow-sample.html"> MeanMean</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>

<p></p>
<p>
Hi everyone, I’m at JSM 2018 right now, so feel free to drop by my session or drop by in the halls! Just give me a <a rel="nofollow" href="https://twitter.com/jlisic" target="_blank">tweet</a>!
</p>
<p>
Back to the meat-and-potatoes of this post. A while ago I was running good old <code>sample</code> and comparing its performance to my <code>lpm2_kdtree</code> function in the <i>BalancedSampling</i> package (Grafström and Lisic, 2016).<br>
In this comparison I noticed that sample was in some cases slower than my balanced sampling method when using sampling weights.
</p>
<div>
<pre><span>library</span><span>(BalancedSampling)</span>
<span>library</span><span>(ggplot2)</span>

<span># sampling proportion</span>
<span>f</span> <span>&lt;-</span> <span>0.1</span>

<span># population size </span>
<span>N_values</span> <span>&lt;-</span>  <span>seq</span><span>(</span><span>from</span><span>=</span><span>5000</span><span>,</span> <span>to</span><span>=</span><span>100000</span><span>,</span><span>length</span><span>=</span><span>10</span><span>)</span>

<span>set.seed</span><span>(</span><span>99</span><span>)</span>

<span>pi</span> <span>&lt;-</span> <span>runif</span><span>(</span><span>100000</span><span>)</span>

<span>result</span> <span>&lt;-</span> <span>c</span><span>()</span>

<span>for</span><span>( N</span> <span>in</span> <span>N_values) {</span>
  <span># get sample size </span>
  <span>n</span> <span>&lt;-</span> <span>round</span><span>(f</span><span>*</span><span>N)</span>

  <span>## run lpm2-kdtree</span>

  <span>start_time</span> <span>&lt;-</span> <span>proc.time</span><span>()[</span><span>3</span><span>]</span>
  <span># to make things even we normalize pi</span>
  <span>pi_tmp</span> <span>&lt;-</span> <span>pi[</span><span>1</span><span>:</span><span>N]</span> <span>*</span> <span>n</span><span>/</span><span>sum</span><span>(pi)</span>
  <span># we probably don't want to balance by pi_tmp, but </span>
  <span># for checking speed this is fine</span>
  <span>lpm2_kdtree</span><span>(pi_tmp,pi_tmp)</span>
  <span>end_time</span> <span>&lt;-</span> <span>proc.time</span><span>()[</span><span>3</span><span>]</span> <span>-</span> <span>start_time</span>

  <span>result</span> <span>&lt;-</span> <span>rbind</span><span>(result,</span> <span>c</span><span>(f, N,</span> <span>1</span><span>, end_time ))</span>

  <span>## run sample </span>

  <span>index</span> <span>&lt;-</span> <span>1</span><span>:</span><span>N</span>

  <span>start_time</span> <span>&lt;-</span> <span>proc.time</span><span>()[</span><span>3</span><span>]</span>
  <span># to make things even we normalize pi</span>
  <span>sample</span><span>(index,</span><span>size</span><span>=n,</span><span>prob</span><span>=pi_tmp)</span>
  <span>end_time</span> <span>&lt;-</span> <span>proc.time</span><span>()[</span><span>3</span><span>]</span> <span>-</span> <span>start_time</span>

  <span>result</span> <span>&lt;-</span> <span>rbind</span><span>(result,</span> <span>c</span><span>(f, N,</span> <span>2</span><span>, end_time ))</span>
<span>}</span>


<span>colnames</span><span>(result)</span> <span>&lt;-</span>
  <span>c</span><span>(</span><span>'SamplingFraction'</span><span>,</span> <span>'PopulationSize'</span><span>,</span> <span>'Method'</span><span>,</span><span>'RunTime'</span><span>)</span>
<span>result</span> <span>&lt;-</span> <span>as.data.frame</span><span>(result)</span>

<span>result[</span><span>c</span><span>(</span><span>'SamplingFraction'</span><span>,</span><span>'Method'</span><span>)]</span> <span>&lt;-</span>
  <span>lapply</span><span>(result[</span><span>c</span><span>(</span><span>'SamplingFraction'</span><span>,</span><span>'Method'</span><span>)], as.factor)</span>
<span>levels</span><span>(result</span><span>$</span><span>Method)</span> <span>&lt;-</span> <span>c</span><span>(</span><span>&quot;lpm2_kdtree&quot;</span><span>,</span><span>&quot;sample&quot;</span><span>)</span>

<span>ggplot</span><span>(</span> <span>aes</span><span>(</span><span>x</span><span>=PopulationSize,</span> <span>y</span><span>=RunTime,</span> <span>class</span><span>=Method),</span> <span>data</span><span>=result)</span><span>+</span>
  <span>geom_line</span><span>(</span><span>aes</span><span>(</span><span>color</span><span>=Method))</span> <span>+</span>
  <span>ylab</span><span>(</span><span>&quot;Run Time (s)&quot;</span><span>)</span> <span>+</span>
  <span>xlab</span><span>(</span><span>&quot;Population Size&quot;</span><span>)</span>
</pre>
</div>
<figure>
<img src="https://i2.wp.com/meanmean.me/blog/sample/figure/unnamed-chunk-2-1.png?w=456" data-recalc-dims="1"><figcaption><b>Figure 1:</b> Lpm2_kdtree vs R sample.</figcaption></figure>
<p> <b>How could this be? </b></p>
<p> Let’s condsider a few possible reasons. </p>
<p> <b>Are you running an old version of R? </b></p>
<p> No, this was done with R 3.5.0, but results hold on 3.5.1.</p>
<p><b> Is this because you are a <i>freedom hating</i> Mac user?</b></p>
<p>Nope, I get similar results in Linux </p>
<p><b>Is this due to the <code>revsort</code> in <code>ProbSampleNoReplace</code> found in <i>random.c</i>?</b></p>
<p>How amazingly prescient of you! If you track down the code path, the method used for sampling has a revsort call right at the begining (<a rel="nofollow" href="https://github.com/wch/r-source/blob/trunk/src/main/random.c" target="_blank">Line 323</a>). Some basic testing seems to identify this as the likely culprit. Therefore, it might be possible to speed up sampling by checking if the probabilities are sorted.</p>
<p>Keep in mind that there are likely good reasons for R to use this method relative to others, including the alternative method below. One is computational stability, split sampling through the pivotal methods tends to work well when sample sizes aren’t extremely small. However, I don’t have any good comparisons of methods, if anyone knows please let me know, doing a literature search on this is pretty far down on my todo list.</p>
<p><b>That’s great, but what can I do to improve my weighted or PPS sampling?</b></p>
<p>Well, I wrote a quick pivotal method implementation (think LPM2 without balancing) that seems to work. It’s available as <code>split_sample</code> on my github version of BalancedSampling.<br>
There is also <code>UPrandompivotal</code> in the <i>sampling</i> package (Tillé and Matei, 2016).<br>
However, my implementation seems to be considerably faster. We will take a look at this in another blog post.
</p>
<p>Finally, keep in mind that split_sample probably isn’t as well tested as <code>UPrandompivotal</code>. So I’m happy to take any feedback. </p>
<div>
<pre><span>library</span><span>(devtools)</span>
<span>install_github</span><span>(</span><span>&quot;jlisic/BalancedSampling&quot;</span><span>)</span>
<span>library</span><span>(BalancedSampling)</span>
<span>library</span><span>(ggplot2)</span>
<span>library</span><span>(sampling)</span>

<span># sampling proportion</span>
<span>f</span> <span>&lt;-</span> <span>0.1</span>

<span># population size </span>
<span>N_values</span> <span>&lt;-</span>  <span>seq</span><span>(</span><span>from</span><span>=</span><span>5000</span><span>,</span> <span>to</span><span>=</span><span>100000</span><span>,</span><span>length</span><span>=</span><span>10</span><span>)</span>

<span>set.seed</span><span>(</span><span>99</span><span>)</span>

<span>pi</span> <span>&lt;-</span> <span>runif</span><span>(</span><span>100000</span><span>)</span>

<span>result</span> <span>&lt;-</span> <span>c</span><span>()</span>

<span>for</span><span>( N</span> <span>in</span> <span>N_values) {</span>
  <span># get sample size </span>
  <span>n</span> <span>&lt;-</span> <span>round</span><span>(f</span><span>*</span><span>N)</span>

  <span>## run lpm2-kdtree</span>

  <span>start_time</span> <span>&lt;-</span> <span>proc.time</span><span>()[</span><span>3</span><span>]</span>
  <span># to make things even we normalize pi</span>
  <span>pi_tmp</span> <span>&lt;-</span> <span>pi[</span><span>1</span><span>:</span><span>N]</span> <span>*</span> <span>n</span><span>/</span><span>sum</span><span>(pi)</span>
  <span># we probably don't want to balance by pi_tmp, but </span>
  <span># for checking speed this is fine</span>
  <span>lpm2_kdtree</span><span>(pi_tmp,pi_tmp)</span>
  <span>end_time</span> <span>&lt;-</span> <span>proc.time</span><span>()[</span><span>3</span><span>]</span> <span>-</span> <span>start_time</span>

  <span>result</span> <span>&lt;-</span> <span>rbind</span><span>(result,</span> <span>c</span><span>(f, N,</span> <span>1</span><span>, end_time ))</span>


  <span>## run pivotal method (BalancedSampling) </span>

  <span>index</span> <span>&lt;-</span> <span>1</span><span>:</span><span>N</span>

  <span>start_time</span> <span>&lt;-</span> <span>proc.time</span><span>()[</span><span>3</span><span>]</span>
  <span># to make things even we normalize pi</span>
  <span>pi_tmp</span> <span>&lt;-</span> <span>pi[</span><span>1</span><span>:</span><span>N]</span> <span>*</span> <span>n</span><span>/</span><span>sum</span><span>(pi)</span>
  <span>split_sample</span><span>(pi_tmp)</span>
  <span>end_time</span> <span>&lt;-</span> <span>proc.time</span><span>()[</span><span>3</span><span>]</span> <span>-</span> <span>start_time</span>

  <span>result</span> <span>&lt;-</span> <span>rbind</span><span>(result,</span> <span>c</span><span>(f, N,</span> <span>2</span><span>, end_time ))</span>


  <span>## run pivotal sampling (sampling) </span>

  <span>index</span> <span>&lt;-</span> <span>1</span><span>:</span><span>N</span>

  <span>start_time</span> <span>&lt;-</span> <span>proc.time</span><span>()[</span><span>3</span><span>]</span>
  <span># to make things even we normalize pi</span>
  <span>pi_tmp</span> <span>&lt;-</span> <span>pi[</span><span>1</span><span>:</span><span>N]</span> <span>*</span> <span>n</span><span>/</span><span>sum</span><span>(pi)</span>
  <span>UPrandompivotal</span><span>(pi_tmp)</span>
  <span>end_time</span> <span>&lt;-</span> <span>proc.time</span><span>()[</span><span>3</span><span>]</span> <span>-</span> <span>start_time</span>

  <span>result</span> <span>&lt;-</span> <span>rbind</span><span>(result,</span> <span>c</span><span>(f, N,</span> <span>3</span><span>, end_time ))</span>
<span>}</span>


<span>colnames</span><span>(result)</span> <span>&lt;-</span>
  <span>c</span><span>(</span><span>'SamplingFraction'</span><span>,</span> <span>'PopulationSize'</span><span>,</span> <span>'Method'</span><span>,</span><span>'RunTime'</span><span>)</span>
<span>result</span> <span>&lt;-</span> <span>as.data.frame</span><span>(result)</span>

<span>result[</span><span>c</span><span>(</span><span>'SamplingFraction'</span><span>,</span><span>'Method'</span><span>)]</span> <span>&lt;-</span>
  <span>lapply</span><span>(result[</span><span>c</span><span>(</span><span>'SamplingFraction'</span><span>,</span><span>'Method'</span><span>)], as.factor)</span>
<span>levels</span><span>(result</span><span>$</span><span>Method)</span> <span>&lt;-</span>
  <span>c</span><span>(</span><span>&quot;lpm2_kdtree&quot;</span><span>,</span><span>&quot;BalancedSampling::split_sample&quot;</span><span>,</span>
    <span>&quot;sampling::UPrandompivotal&quot;</span><span>)</span>

<span>ggplot</span><span>(</span> <span>aes</span><span>(</span><span>x</span><span>=PopulationSize,</span> <span>y</span><span>=RunTime,</span> <span>class</span><span>=Method),</span> <span>data</span><span>=result)</span><span>+</span>
  <span>geom_line</span><span>(</span><span>aes</span><span>(</span><span>color</span><span>=Method))</span> <span>+</span>
  <span>ylab</span><span>(</span><span>&quot;Run Time (s)&quot;</span><span>)</span> <span>+</span>
  <span>xlab</span><span>(</span><span>&quot;Population Size&quot;</span><span>)</span>
</pre>
</div>
<figure>
<img src="https://i2.wp.com/meanmean.me/blog/sample/figure/unnamed-chunk-3-1.png?w=456" data-recalc-dims="1"><figcaption><b>Figure 2:</b> Lpm2_kdtree vs BalancedSampling::split_sample vs sampling::UPrandompivotal.</figcaption></figure>
<p> Have a great JSM! </p>
<h2>References</h2>
<p>
Grafström, A., Lisic, J.: BalancedSampling: balanced and spatially balanced sampling. R package version 1.5.2. <u>https://CRAN.R-project.org/package=BalancedSampling</u>. Accessed August 28, 2018. (2016)
</p>
<p>
Tillé, Y., Matei, M.: sampling: Survey Sampling. R package version 2.8. <u>https://CRAN.R-project.org/package=sampling</u>. Accessed August 28, 2018. (2016)</p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://meanmean.me/2018/07/28/Mysteriously-slow-sample.html"> MeanMean</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/X_eBmUuRfOE" height="1"></a></div>