Reputation: 557
I have a question on manipulating data in data.frame.
Essentially I have a large data set - abbreviated version below:
structure(list(nm_mean = c(194213914.326, 194213914.326, 194213914.326,
194213914.326, 194213914.326, 217947112.739), nm_se = c(9984735.05918367,
9984735.05918367, 9984735.05918367, 9984735.05918367, 9984735.05918367,
11010386.0760204), alpha = c(193.197697846336, 214.592588477741,
240.246557258741, 258.116959355425, 282.560024775668, 306.610038660465
), beta = c(61526.2664158025, 57950.9563448233, 56085.1512614369,
52919.4794239927, 51483.4591654126, 50405.8186695088)), .Names = c("nm_mean",
"nm_se", "alpha", "beta"), row.names = c(NA, 6L), class = "data.frame")
I want to use rbeta to generate probabilities using the beta distribution and alpha and beta as the parameters
Similarly I want to use rnorm to generate random numbers using the normal distribution with nm_mean and nm_se as the mean and sd.
I then want to multiply the rbeta values generated by the rnorm values and extract the 50th, 25th and 75th quantile back into the dataframe
So as an example for row 1
x <- rbeta(1000,193.1977,61526.27)
y <- rnorm(1000,194213914,9984735)
z <- x*y
dat$ce <- quantile(z,0.5)
dat$ll <- quantile(z,0.25)
dat$ul <- quantile(z,0.975)
In essence i get a ce, ll and ul for product of the rbeta and rnorm appended back to the database.
Upvotes: 3
Views: 224
Reputation: 23241
This is vectorized solution based on my conversation with @thelatemail:
n <- 1000
grp <- nrow(dat)
z <- with(dat, rnorm(grp*n, nm_mean, nm_se) * rbeta(grp*n, alpha, beta) )
m <- 1
for(i in 1:nrow(dat)){
dat$ce[i] <- quantile(z[m:(i*1000)],0.5)
dat$ll[i] <- quantile(z[m:(i*1000)],0.25)
dat$ul[i] <- quantile(z[m:(i*1000)],0.975)
m <- m + 1000
}
A less vectorized solution is:
for(i in 1:nrow(dat)){
x <- rbeta(1000, shape1 = dat$alpha[i], shape2 = dat$beta[i])
y <- rnorm(n=1000,dat$nm_mean[i],dat$nm_se[i])
z <- x*y
dat$ce[i] <- quantile(z,0.5)
dat$ll[i] <- quantile(z,0.25)
dat$ul[i] <- quantile(z,0.975)
}
dat
nm_mean nm_se alpha beta ce ll ul 1 194213914 9984735 193.1977 61526.27 607563.9 573229.9 713057.2 2 194213914 9984735 214.5926 57950.96 712268.5 674826.3 836950.8 3 194213914 9984735 240.2466 56085.15 823322.9 777482.8 981156.7 4 194213914 9984735 258.1170 52919.48 937331.2 884945.0 1095876.3 5 194213914 9984735 282.5600 51483.46 1059980.4 1003596.4 1225615.6 6 217947113 11010386 306.6100 50405.82 1316733.1 1250190.1 1515185.0
Upvotes: 1
Reputation: 161155
Motivated by @HackR's code, what I think is a functional vectorized version:
set.seed(42)
n <- 1000
nrows <- nrow(dat)
rn <- matrix(rnorm(nrows * n, dat$nm_mean, dat$nm_se), ncol = nrows, byrow = TRUE)
rb <- matrix(rbeta(nrows * n, shape1 = dat$alpha, shape2 = dat$beta),
ncol = nrows, byrow = TRUE)
cbind(dat,
structure(t(apply(rn * rb, 2, function(z) quantile(z, c(0.5, 0.25, 0.975)))),
.Dimnames = list(NULL, c("ce", "ll", "ul"))))
# nm_mean nm_se alpha beta ce ll ul
# 1 194213914 9984735 193.1977 61526.27 608455.3 570100.5 710373.6
# 2 194213914 9984735 214.5926 57950.96 715305.0 677754.3 856570.7
# 3 194213914 9984735 240.2466 56085.15 825143.7 778351.2 979361.1
# 4 194213914 9984735 258.1170 52919.48 943261.4 895832.6 1091899.3
# 5 194213914 9984735 282.5600 51483.46 1054514.3 995640.8 1226176.4
# 6 217947113 11010386 306.6100 50405.82 1312325.0 1247030.8 1515630.5
Upvotes: 6