[R] problem with outer

[Joshua Wiley]

Great, that helps.  With all those variables, your function worked for
me.  I did slightly edit it for testing, mostly just adding spaces to
make it easier for me to read.


N <- 10; q <- 8; k <- q-2; L <- 5; n <- 2 * k + L + 1
p_0 <- c(0.2, 0.4, 0.4)
f_1 <- qf(0.05, 2, 2)
f_2 <- qf(0.95, 2, 2)
p_11 <- seq(0, 0.3, 0.1)
p_12 <- seq(0.1, 0.4, 0.1)

guete <- function(p_11, p_12) {
  set.seed(1000)
  S_vek <- matrix(0, nrow = N, ncol = 1)
  zz <- cbind(p_11, p_12, (1 - p_11 - p_12))

  for(i in 1:N) {
    X_0 <- rmultinom(q - 1, size = 1, prob = p_0)
    X_1 <- rmultinom(n - q + 1, size = 1, prob = zz)
    N_0 <- apply(X_0[ ,(n - 2 * k - L + 1):(n - k - L)], 1, sum)
    N_1 <- apply(X_1[ ,(n - q - k + 2):(n - q + 1)], 1, sum)
    S_vek[i] <- ( (sum( ((N_1 - k * zz)^2) / k * zz)) /
                 (sum( ((N_0 - k * p_0)^2) / k * p_0))
                 ) - 1
  }

  1 - mean(f_1 <= S_vek & S_vek <= f_2)
}

guete(p_11 = p_11, p_12 = p_12)


For your second problem:

p_11=seq(0,1,0.1)
p_12=seq(0,1,0.1)
then i get also this error message:
Error in rmultinom(n - q + 1, size = 1, prob = rbind(p_11, p_12, (1 -  :
 non-positive probability


the probabilities for rmultinom need to be non-negative.  But when you
do seq(0, 1, .1), you get some 1s in p_11 and p_12.  This means part
of your definition for the probs argument is (1 - 1 - 1) = -1.  So you
need to make sure that that maximum value of p_11 + p_12 is <= 1 (if
you are going to subtract both from 1 and use that as a probability).

Hope that helps,

Josh



On Wed, Sep 8, 2010 at 3:24 PM, tuggi <tugbagueclue at web.de> wrote:
>
> hello,
>
> thank you for your answer. here the parameter
>
>  N=10
>
>
>  q=8
>
>
>
>  k=q-2
>
>
>
>  L=5
>
>
>
>  n=2*k+L+1
>
>
>  p_0=c(0.2,0.4,0.4)
>
> f_1=qf(0.05,2,2)
>
>
> f_2=qf(0.95,2,2)
>
> tuggi
> --
> View this message in context: http://r.789695.n4.nabble.com/problem-with-outer-tp2532074p2532132.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.com/