title: How to create a more efficient simulation loop for Monte Carlo in R
body: <p>The purpose of this exercise is to create a population distribution of nutrient intake values. There were repeated measures in the earlier data, these have been removed so each row is a unique person in the data frame.</p>
  
  <p>I have this code, which works quite well when tested with a small number of my data frame rows. For all 7135 rows, it is very slow. I tried to time it, but I crashed it out when the elapsed running time on my machine was 15 hours. The <code>system.time</code> results were <code>Timing stopped at: 55625.08 2985.39 58673.87</code>.</p>
  
  <p>I would appreciate any comments on speeding up the simulation:</p>
  
  <pre><code>Male.MC &lt;-c()
  for (j in 1:100)            {
  for (i in 1:nrow(Male.Distrib))  {
      u2        &lt;- Male.Distrib$stddev_u2[i] * rnorm(1, mean = 0, sd = 1)
      mc_bca    &lt;- Male.Distrib$FixedEff[i] + u2
      temp      &lt;- Lambda.Value*mc_bca+1
      ginv_a    &lt;- temp^(1/Lambda.Value)
      d2ginv_a  &lt;- max(0,(1-Lambda.Value)*temp^(1/Lambda.Value-2))
      mc_amount &lt;- ginv_a + d2ginv_a * Male.Resid.Var / 2
  z &lt;- data.frame(
       RespondentID = Male.Distrib$RespondentID[i], 
       Subgroup     = Male.Distrib$Subgroup[i], 
       mc_amount    = mc_amount,
       IndvWeight   = Male.Distrib$INDWTS[i]/100
       )
  
  Male.MC &lt;- as.data.frame(rbind(Male.MC,z))
      }
  }
  </code></pre>
  
  <p>For each of the 7135 observations in my dataset, 100 simulated nutrient values are created, then back transformed to the original measurement level (the simulation is using the results from a nonlinear mixed effect model on BoxCox transformed nutrient values).</p>
  
  <p>I would prefer not to use <code>for</code> loops, as I read that they are inefficient in <code>R</code> but I do not understand enough about options based on <code>apply</code> to use those as an alternative. <code>R</code> is being run on stand-alone machines, normally this would be a standard Dell-type desktop running a Windows 7 variant, if that influences the recommendations for how to change the code.</p>
  
  <p>Update: To reproduce this for testing,
  <code>Lambda.Value</code>=0.4 and <code>Male.Resid.Var</code>=12.1029420429778 and <code>Male.Distrib$stddev_u2</code> is a constant value over all observations.</p>
  
  <p><code>str(Male.Distrib)</code> is</p>
  
  <pre><code>'data.frame':   7135 obs. of  14 variables:
   $ RndmEff     : num  1.34 -5.86 -3.65 2.7 3.53 ...
   $ RespondentID: num  9966 9967 9970 9972 9974 ...
   $ Subgroup    : Ord.factor w/ 6 levels "3"&lt;"4"&lt;"5"&lt;"6"&lt;..: 4 3 2 4 1 4 2 5 1 2 ...
   $ RespondentID: int  9966 9967 9970 9972 9974 9976 9978 9979 9982 9993 ...
   $ Replicates  : num  41067 2322 17434 21723 375 ...
   $ IntakeAmt   : num  33.45 2.53 9.58 43.34 55.66 ...
   $ RACE        : int  2 3 2 2 3 2 2 2 2 1 ...
   $ INDWTS      : num  41067 2322 17434 21723 375 ...
   $ TOTWTS      : num  1.21e+08 1.21e+08 1.21e+08 1.21e+08 1.21e+08 ...
   $ GRPWTS      : num  41657878 22715139 10520535 41657878 10791729 ...
   $ NUMSUBJECTS : int  1466 1100 1424 1466 1061 1466 1424 1252 1061 1424 ...
   $ TOTSUBJECTS : int  7135 7135 7135 7135 7135 7135 7135 7135 7135 7135 ...
   $ FixedEff    : num  6.09 6.76 7.08 6.09 6.18 ...
   $ stddev_u2   : num  2.65 2.65 2.65 2.65 2.65 ...
  </code></pre>
  
  <p><code>head(Male.Distrib)</code> is</p>
  
  <pre><code>    RndmEff RespondentID Subgroup RespondentID Replicates IntakeAmt RACE INDWTS    TOTWTS   GRPWTS NUMSUBJECTS TOTSUBJECTS  FixedEff stddev_u2
  1  1.343753         9966        6         9966      41067 33.449808    2  41067 120622201 41657878        1466        7135  6.089918  2.645938
  2 -5.856516         9967        5         9967       2322  2.533528    3   2322 120622201 22715139        1100        7135  6.755664  2.645938
  3 -3.648339         9970        4         9970      17434  9.575439    2  17434 120622201 10520535        1424        7135  7.079757  2.645938
  4  2.697533         9972        6         9972      21723 43.340180    2  21723 120622201 41657878        1466        7135  6.089918  2.645938
  5  3.531878         9974        3         9974        375 55.660607    3    375 120622201 10791729        1061        7135  6.176319  2.645938
  6  6.627767         9976        6         9976      48889 91.480049    2  48889 120622201 41657878        1466        7135  6.089918  2.645938
  </code></pre>
  
  <p>Update 2: the line of the function that is causing the <code>NaN</code> results is </p>
  
  <pre><code>d2ginv_a  &lt;- max(0,(1-Lambda.Value)*temp^(1/Lambda.Value-2))
  </code></pre>
  
  <p>Thanks to everyone for their assistance and comments, and also for the speed of responses.</p>
  
  <p>Update: @Ben Bolker is correct that it is the negative <code>temp</code> values that are causing the NaN issue. I missed this with some testing (after commenting out the function so that only the <code>temp</code> values are returned, and calling my result data frame <code>Test</code>). This code reproduces the <code>NaN</code> issue:</p>
  
  <pre><code>&gt; min(Test)
  [1] -2.103819
  &gt; min(Test)^(1/Lambda.Value)
  [1] NaN
  </code></pre>
  
  <p>But putting the value in as a value and then running the same(?) calculation gives me a result, so I missed this when doing manual calculations:</p>
  
  <pre><code>&gt; -2.103819^(1/Lambda.Value) 
  [1] -6.419792
  </code></pre>
  
  <p>I now have working code that (I think) uses vectorization, and it is blindingly fast. Just in case anyone else has this issue, I am posting the working code below. I've had to add a minimum to prevent the &lt;0 issue with the calculation. Thank you to everyone who helped, and to coffee. I did try putting the <code>rnorm</code> results to a dataframe, and that really slowed things down, creating them this way and then using <code>cbind</code> is really quick. <code>Male.Distrib</code> is my full data frame of 7135 observations, but this code should work on the cutdown version I posted earlier (not tested).</p>
  
  <pre><code>Min_bca &lt;- ((.5*min(Male.AddSugar$IntakeAmt))^Lambda.Value-1)/Lambda.Value
  Test &lt;- Male.Distrib[rep(seq.int(1,nrow(Male.Distrib)), 100), 1:ncol(Male.Distrib)]
  RnormOutput &lt;- rnorm(nrow(Test),0,1)
  Male.Final &lt;- cbind(Test,RnormOutput)
  Male.Final$mc_bca    &lt;- Male.Final$FixedEff + (Male.Final$stddev_u2 *     Male.Final$RnormOutput)
  Male.Final$temp      &lt;- ifelse(Lambda.Value*Male.Final$mc_bca+1 &gt; Lambda.Value*Min_bca+1,
                             Lambda.Value*Male.Final$mc_bca+1, Lambda.Value*Min_bca+1)
  Male.Final$ginv_a    &lt;- Male.Final$temp^(1/Lambda.Value)
  Male.Final$d2ginv_a  &lt;- ifelse(0 &gt; (1-Lambda.Value)*Male.Final$temp^(1/Lambda.Value-2),
                             0, (1-Lambda.Value)*Male.Final$temp^(1/Lambda.Value-2))
  Male.Final$mc_amount &lt;- Male.Final$ginv_a + Male.Final$d2ginv_a * Male.Resid.Var / 2
  </code></pre>
  
  <p>Lessons for the day:</p>
  
  <ul>
  <li>a distribution function does not appear to be resampled in a loop if you try to do what I was trying earlier</li>
  <li>you can't use <code>max()</code> the way I tried, as it returns the maximum value from the column, whereas I wanted the maximum from two values. The <code>ifelse</code> statement is the replacement one to do.</li>
  </ul>