The defintion of the sequency of probability was wrong!
> p_11=seq(0,1,0.1)
> p_12=seq(0,1,0.1)
Since your multinomial distribution requires P_11, P_12, P13=1-P_11-P12 be
greater than or equal to zero and P_11+P_12+P_13 = 1. In your above
definition of P_11[6] =P_12[6]= 0.6, P_11[7] = P_12[7] = 0.7, ..., P_11[10]
=P_12[10]=1.0. For these pairs 1-p_11-p_12 < 0! There is why you had the
error message!
> X_1=rmultinom(n-q+1,size=1,prob=cbind(p_11,p_12,(1-p_11-p_12)))
Suggested Solution:
Redefine P_11 and P_12 such that P_11[i] + P_12[i] <1 for all i! This will
sove your problem!
# you choose whatever a and b such that a+b<=1. for example,
a = 0.4
b = 0.6
# you need to adjust the increaments to get vectors with equal length.
p_11=seq(0,a,0.05*a)
p_12=seq(0, b,0.05*b)
>
>
> guete = function(p_11,p_12) {
+ if(p_11+p_12<1)
+ set.seed(1000)
+ S_vek=matrix(0,nrow=N,ncol=1)
+ for(i in 1:N) {
+ X_0=rmultinom(q-1,size=1,prob=p_0)
+
X_1=rmultinom(n-q+1,size=1,prob=cbind(p_11,p_12,(1-p_11-p_12)))
+ 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*cbind(p_11,p_12,(1-p_11-p_12)))^2)/k*cbind(p_11,p_12,(1-p_11-p_12))))/(sum(((N_0-k*p_0)^2)/k*p_0)))-1
+ }
+ 1-mean(f_1<=S_vek & S_vek <=f_2)
+ }
>
>
> f=outer(p_11,p_12,Vectorize(guete))
> f
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
[,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
[1,] 1.0 0.9 0.9 0.8 0.8 0.7 0.6 0.6 0.6 0.4 0.7 0.8 0.8
0.7 0.6 0.7 0.9 0.8 0.7 0.5 0.6
[2,] 0.9 1.0 0.9 0.7 0.6 0.6 0.5 0.5 0.7 0.7 0.6 0.5 0.4
0.2 0.3 0.4 0.7 0.6 0.2 0.3 0.3
[3,] 0.9 1.0 0.8 0.8 0.7 0.7 0.7 0.6 0.6 0.5 0.7 0.6 0.3
0.3 0.4 0.5 0.5 0.5 0.4 0.3 0.3
[4,] 0.9 0.8 0.8 0.8 0.7 0.7 0.5 0.5 0.6 0.8 0.7 0.6 0.5
0.5 0.4 0.6 0.5 0.6 0.5 0.4 0.5
[5,] 0.8 0.8 0.9 0.7 0.6 0.5 0.4 0.4 0.5 0.6 0.4 0.6 0.5
0.5 0.5 0.5 0.6 0.4 0.4 0.6 0.5
[6,] 0.8 0.9 0.7 0.5 0.6 0.6 0.5 0.4 0.5 0.6 0.5 0.6 0.3
0.4 0.6 0.5 0.6 0.4 0.3 0.6 0.5
[7,] 0.8 1.0 0.6 0.4 0.5 0.7 0.6 0.6 0.7 0.7 0.5 0.4 0.4
0.4 0.5 0.6 0.4 0.4 0.4 0.6 0.6
[8,] 0.7 1.0 0.6 0.5 0.5 0.7 0.7 0.8 0.8 0.7 0.6 0.4 0.6
0.6 0.7 0.7 0.6 0.5 0.5 0.5 0.8
[9,] 0.7 0.7 0.6 0.7 0.5 0.4 0.5 0.6 0.6 0.4 0.2 0.2 0.5
0.5 0.5 0.5 0.5 0.2 0.2 0.4 0.6
[10,] 0.6 0.7 0.7 0.5 0.4 0.4 0.5 0.7 0.7 0.5 0.3 0.4 0.5
0.4 0.5 0.4 0.5 0.3 0.3 0.6 0.7
[11,] 0.7 0.7 0.8 0.7 0.6 0.5 0.5 0.7 0.7 0.5 0.5 0.6 0.5
0.5 0.5 0.5 0.5 0.5 0.5 0.7 0.7
[12,] 0.7 0.7 0.6 0.5 0.5 0.4 0.5 0.5 0.6 0.5 0.4 0.2 0.2
0.3 0.2 0.2 0.4 0.5 0.6 0.5 0.5
[13,] 0.6 0.6 0.6 0.4 0.3 0.2 0.7 0.6 0.5 0.6 0.4 0.4 0.3
0.3 0.4 0.4 0.6 0.6 0.6 0.7 0.3
[14,] 0.4 0.6 0.6 0.4 0.1 0.7 0.7 0.5 0.6 0.6 0.5 0.4 0.3
0.3 0.4 0.5 0.6 0.5 0.6 0.6 0.2
[15,] 0.4 0.6 0.6 0.4 0.3 0.4 0.4 0.3 0.4 0.5 0.5 0.5 0.3
0.5 0.5 0.5 0.4 0.3 0.4 0.4 0.3
[16,] 0.7 0.3 0.3 0.2 0.2 0.2 0.1 0.1 0.2 0.1 0.1 0.2 0.2
0.1 0.1 0.2 0.1 0.1 0.0 0.3 0.2
[17,] 0.7 0.4 0.4 0.2 0.2 0.2 0.2 0.2 0.2 0.1 0.3 0.3 0.3
0.3 0.1 0.2 0.2 0.2 0.2 0.2 0.2
[18,] 0.8 0.3 0.3 0.2 0.2 0.2 0.1 0.1 0.1 0.1 0.1 0.1 0.1
0.1 0.1 0.1 0.1 0.2 0.2 0.2 0.3
[19,] 0.8 0.7 0.7 0.6 0.5 0.4 0.4 0.3 0.2 0.2 0.2 0.2 0.2
0.1 0.3 0.4 0.4 0.5 0.5 0.7 0.7
[20,] 0.6 0.4 0.3 0.2 0.3 0.3 0.3 0.3 0.2 0.3 0.3 0.3 0.2
0.3 0.3 0.3 0.3 0.3 0.2 0.3 0.4
[21,] 0.6 0.3 0.4 0.4 0.4 0.3 0.3 0.4 0.4 0.1 0.1 0.1 0.4
0.4 0.3 0.3 0.4 0.4 0.4 0.3 0.5
--
View this message in context:
http://r.789695.n4.nabble.com/problem-with-outer-tp2532074p2532282.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
R-help@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.
