Saving quantities is a for loop

Question

I am having problems when saving the results in a for loop. I am computing a variance (this is not relevant I think) and my code is:

library(dirmult)
n <- 50
p <- 20
size <- 5*p
prob_true <- rep(1/p, p)
multinom <- as.matrix(rmultinom(n, size, prob = prob_true))
zeros <- round(0.5*p*n)
a <- c(as.matrix(multinom))
a[sample(1:(p*n), zeros)] <- 0
data_zeros <- matrix(a, p, n)
dirmult <- dirmult(t(data_zeros))
alpha <- dirmult$gamma
sum_alpha <- (1-dirmult$theta)/dirmult$theta

for (j in ncol(data_zeros)){
  A <- alpha/sum_alpha
  B <- 1 - A
  N <- colSums(data_zeros)
  C <- 1 + sum_alpha
  var_s_dirm <- list()
  var_s_dirm[[j]] <- N[j]*A*B*((N[j]+sum_alpha)/C)
}

In particular I can say that alpha is a vector with 20 values, sum_alpha is a scalar data_zeros is my dataset which has 20 rows and 50 columns and N is the sum of each column of the dataset, so it is a vector with 50 values.

It seems very simple to do what I wanted to do: I want to get a list with 50 vectors where each one differs form the other by the fact that I multiply for a different value of N.

I really hope that somebody can help me finding the error.

Hard to find/fix errors if we cannot reproduce the problem.. Please add sample data.. — Wimpel, Mar 04 '21 at 18:58
Your code till does not run.. `Error in rmultinom(n, size, prob = prob_true) : object 'size' not found`. it's probably a good idea to read how to provide an minimal reproducible example: https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example — Wimpel, Mar 04 '21 at 19:10
Yes, sorry I forgot that line of code, my bad! I have many lines and I got confused — Bibi, Mar 04 '21 at 19:12
That is probably why it is best practive to use the guidelines in the link above ;-) — Wimpel, Mar 04 '21 at 19:13
You're right! For now on I will look at it before asking questions... but now I think that this one works fine... — Bibi, Mar 04 '21 at 19:15
still does not run.. `Error in dirmult(t(data_zeros[[i]])) : could not find function "dirmult"` please try your code in a clean R session before you post it. Also: you will probably get faster/better answers if you use the guides in the link above. — Wimpel, Mar 04 '21 at 19:15
also, from *the R inferno*, chapter 9: It can not be emphasized enough that a message is putting yourself at the mercy of strangers. If someone has the wit and knowledge to answer your question, they probably have other things they would like to do. Making your message clear, concise and user-friendly gives you the best hope of at least one of those strangers diverting their attention away from their life towards your problem. — Wimpel, Mar 04 '21 at 19:18
Done! Now I get in var_s_dirm all NA except the last quantity... that's the problem.. — Bibi, Mar 04 '21 at 19:19
Yes, you're right, I am very sorry... I am new to everything and also coding... this is why sometimes I am not precise, I do no think at all the things that are necessary... sorry again! — Bibi, Mar 04 '21 at 19:24

score 1 · Accepted Answer · answered Mar 04 '21 at 19:29

The problem is (probably) you are setting constants in each time j is increased, and in each step you clear the list with the line var_s_dirm <- list()...

See if this works for you

library(dirmult)
n <- 50
p <- 20
size <- 5*p
prob_true <- rep(1/p, p)
multinom <- as.matrix(rmultinom(n, size, prob = prob_true))
zeros <- round(0.5*p*n)
a <- c(as.matrix(multinom))
a[sample(1:(p*n), zeros)] <- 0
data_zeros <- matrix(a, p, n)
dirmult <- dirmult(t(data_zeros))
alpha <- dirmult$gamma
sum_alpha <- (1-dirmult$theta)/dirmult$theta


A <- alpha/sum_alpha
B <- 1 - A
N <- colSums(data_zeros)
C <- 1 + sum_alpha
var_s_dirm <- list()

for (j in 1:ncol(data_zeros)){
  var_s_dirm[[j]] <- N[j]*A*B*((N[j]+sum_alpha)/C)
}

output

var_s_dirm
[[1]]
 [1] 2.614833 2.327105 2.500483 3.047700 2.233528 2.130223 2.700103 2.869699 2.930213 2.575903 2.198459 2.846096
[13] 2.425448 3.517559 3.136266 2.565345 2.578267 2.763113 2.709707 3.420792

[[2]]
 [1] 2.568959 2.286279 2.456615 2.994231 2.194343 2.092850 2.652732 2.819353 2.878806 2.530712 2.159889 2.796165
[13] 2.382897 3.455848 3.081244 2.520339 2.533034 2.714637 2.662168 3.360778

[[3]]
 [1] 3.211199 2.857849 3.070769 3.742790 2.742930 2.616064 3.315916 3.524193 3.598509 3.163391 2.699862 3.495207
[13] 2.978622 4.319811 3.851556 3.150424 3.166294 3.393297 3.327711 4.200974

....

Saving quantities is a for loop

1 Answers1