StackOverflow2013

Note that there are some explanatory texts on larger screens.

plurals

PO
text
Body
copied!<p>Here is an approach that addresses the 2 biggest speed issues:</p> <ol> <li>Instead of looping over observations(<code>i</code>), we compute them all at once.</li> <li>Instead of looping over MC replications (<code>j</code>), we use <code>replicate</code>, which is a simplified <code>apply</code> meant for this purpose.</li> </ol> <p>First we load the dataset and define a function for what you were doing.</p> <pre><code>Male.Distrib = read.table('MaleDistrib.txt', check.names=F) getMC <- function(df, Lambda.Value=0.4, Male.Resid.Var=12.1029420429778) { u2 <- df$stddev_u2 * rnorm(nrow(df), mean = 0, sd = 1) mc_bca <- df$FixedEff + u2 temp <- Lambda.Value*mc_bca+1 ginv_a <- temp^(1/Lambda.Value) d2ginv_a <- max(0,(1-Lambda.Value)*temp^(1/Lambda.Value-2)) mc_amount <- ginv_a + d2ginv_a * Male.Resid.Var / 2 mc_amount } </code></pre> <p>Then we replicate it a bunch of times.</p> <pre><code>> replicate(10, getMC(Male.Distrib)) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 36.72374 44.491777 55.19637 23.53442 23.260609 49.56022 31.90657 25.26383 25.31197 20.58857 [2,] 29.56115 18.593496 57.84550 22.01581 22.906528 22.15470 29.38923 51.38825 13.45865 21.47531 [3,] 61.27075 10.140378 75.64172 28.10286 9.652907 49.25729 23.82104 31.77349 16.24840 78.02267 [4,] 49.42798 22.326136 33.87446 14.00084 25.107143 25.75241 30.20490 33.14770 62.86563 27.33652 [5,] 53.45546 9.673162 22.66676 38.76392 30.786100 23.42267 28.40211 35.95015 43.75506 58.83676 [6,] 34.72440 23.786004 63.57919 8.08238 12.636745 34.11844 14.88339 21.93766 44.53451 51.12331 </code></pre> <p>Then you can reformat, add IDs, etc., but this is the idea for the main computational part. Good luck!</p>

Querying!

Guidance

An individual column

Larger individual text columns get their own page to allow for proper reading.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload