[R] problem with outer

Peng, C cpeng.usm at gmail.com
Thu Sep 9 04:21:07 CEST 2010


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.



More information about the R-help mailing list