[R] max row

arun smartpink111 at yahoo.com
Sun Mar 10 04:53:17 CET 2013



HI,

Using
c11<- 0.01
c12<- 0.01
c1<- 0.10
c2<- 0.10


One possible problem is that:
dim(res5)
#[1] 513  20

res6<-aggregate(.~m1+n1+m+n,data=res5[,c(1:6,9:12,21:24)] ,max) 
#Error in `[.data.frame`(res5, , c(1:6, 9:12, 21:24)) : 
 # undefined columns selected
A.K.




________________________________
From: Joanna Zhang <zjoanna2013 at gmail.com>
To: arun <smartpink111 at yahoo.com> 
Sent: Saturday, March 9, 2013 10:38 PM
Subject: Re: max row


I think something is wrong with EN, it should not be integer. Let me check.



On Sat, Mar 9, 2013 at 9:36 PM, arun <smartpink111 at yahoo.com> wrote:

In the example set you gave, both N and EN are the same. that means, there are multiple miniumn values for both N and EN.
>
>
>
>
>
>
>
>________________________________
>From: Joanna Zhang <zjoanna2013 at gmail.com>
>To: arun <smartpink111 at yahoo.com>
>Sent: Saturday, March 9, 2013 10:34 PM
>
>Subject: Re: max row
>
>
>if multiple rows with minimum N, I want to get the rows with minimum value of EN.
>
>
>
>On Sat, Mar 9, 2013 at 9:32 PM, arun <smartpink111 at yahoo.com> wrote:
>
>I hope your problem is solved.
>>
>>
>>
>>
>>
>>
>>
>>________________________________
>>From: Joanna Zhang <zjoanna2013 at gmail.com>
>>To: arun <smartpink111 at yahoo.com>
>>Sent: Saturday, March 9, 2013 10:20 PM
>>
>>Subject: Re: max row
>>
>>
>>I know, i don't know how to use the seq in lapply, I tried to use this, but want to try lapply, it may more efficient.
>>
>>
>>for (c11 in seq(0.05,0.06, by=0.01)) {
>>for (c12 in seq(0.05,0.06, by=0.01)) {
>>for (c1 in seq(0.2,0.3,by=0.1)){
>>for (c2 in seq(0.2,0.3,by=0.1)){
>>search(c11,c12,c1,c2)
>>}}}}
>>
>>
>>
>>
>>
>>On Sat, Mar 9, 2013 at 9:17 PM, arun <smartpink111 at yahoo.com> wrote:
>>
>>Hi,
>>>This is deadset wrong.
>>>
>>>Here c11, c12, are only taking the value 0.01
>>>
>>> lapply(0.01:0.03,function(c11) c11)
>>>[[1]]
>>>[1] 0.01
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>________________________________
>>>From: Joanna Zhang <zjoanna2013 at gmail.com>
>>>To: arun <smartpink111 at yahoo.com>
>>>Sent: Saturday, March 9, 2013 10:04 PM
>>>
>>>Subject: Re: max row
>>>
>>>
>>>great. I want to find the rows with the minimum value of EN in one dataset fopt) and the rows with minimum of N in another data set(fmin) , and then I need to loop through different combinations of c11,c12,c1,c2,, c11, c12 need to take values of 0.01,0.02.0.03, and c1, c2 take vaules of 0.10, 0.20, 0.30. And then find the rows with minimum value of EN and the rows with minimum value of N from all the fopt and fmin 
>>>
>>>
>>>final<-do.call(rbind,lapply(0.01:0.03,function(c11)
>>>do.call(rbind,lapply(0.01:0.03,function(c12)
>>>do.call(rbind,lapply(0.10:0.30,function(c1)
>>>do.call(rbind,lapply(0.10:0.30,function(c2)
>>>search (c11,c12,c1,c2) )) )) )) ))
>>>final
>>>
>>>
>>>
>>>
>>>On Sat, Mar 9, 2013 at 8:57 PM, arun <smartpink111 at yahoo.com> wrote:
>>>
>>>After running this:
>>>>I am getting:
>>>>
>>>>
>>>>search(0.01,0.02,0.15,0.20)
>>>>
>>>>  m1 n1 m n cterm2_P0L cterm2_P0H N EN         BH        BL        AH        AL
>>>>1  2  2 4 4  0.8145062  0.6634204 8  8 0.10011292 0.3164062 0.3365796 0.1854938
>>>>2  2  2 5 4  0.7737809  0.7961045 9  9 0.20022583 0.2373047 0.2038955 0.2262191
>>>>3  3  2 5 4  0.7737809  0.7961045 9  9 0.20022583 0.2373047 0.2038955 0.2262191
>>>>4  2  2 4 5  0.8145062  0.6302494 9  9 0.07508469 0.3164063 0.3697506 0.1854938
>>>>5  2  3 4 5  0.8145062  0.6302494 9  9 0.07508469 0.3164063 0.3697506 0.1854938
>>>>
>>>>
>>>>So, here, what you want to do? To find the minimum value of `AL`? 
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>________________________________
>>>>From: Joanna Zhang <zjoanna2013 at gmail.com>
>>>>To: arun <smartpink111 at yahoo.com>
>>>>Sent: Saturday, March 9, 2013 9:49 PM
>>>>
>>>>Subject: Re: max row
>>>>
>>>>
>>>>ok. this is the codes that I used for small data and without the parameter loop, but would have CPU problem if I use the real data. Please run it and let me know if it works for you.
>>>>
>>>>search<-function(c11,c12,c1,c2){
>>>>d<-do.call(rbind,lapply(2:(maxN-6),function(m1) 
>>>>do.call(rbind,lapply(2:(maxN-m1-4),function(n1)
>>>>do.call(rbind,lapply(0:m1,function(x1)
>>>>do.call(rbind,lapply(0:n1,function(y1)
>>>>expand.grid(m1,n1,x1,y1))))))))) 
>>>>names(d)<-c("m1","n1","x1","y1")
>>>>
>>>>d<-within(d,{
>>>>term1_p1<- dbinom(x1,m1, p1L, log=FALSE)* dbinom(y1,n1,p1H, log=FALSE)
>>>>term1_p0<- dbinom(x1,m1, p0L, log=FALSE)* dbinom(y1,n1,p0H, log=FALSE)
>>>>})      
>>>>
>>>>########## add Qm Qn ##################################
>>>>set.seed(8)
>>>>d1<-do.call(rbind,lapply(seq_len(nrow(d)),function(i){
>>>>Pm<- rbeta(1000,0.2+d[i,"x1"],0.8+d[i,"m1"]-d[i,"x1"]);
>>>>Pn<- rbeta(1000,0.2+d[i,"y1"],0.8+d[i,"n1"]-d[i,"y1"]);
>>>>Fm<- ecdf(Pm);
>>>>Fn<- ecdf(Pn);
>>>>
>>>>Fmm<- Fm(p1L);
>>>>Fnn<- Fn(p1H);
>>>>R<- (Fmm+Fnn)/2; 
>>>>Fmm_f<- max(R, Fmm);
>>>>Fnn_f<- min(R, Fnn);
>>>>Qm<- 1-Fmm_f;
>>>>Qn<- 1-Fnn_f;
>>>>data.frame(Qm,Qn)}))
>>>>d2<-cbind(d,d1)
>>>>head(d2)
>>>>
>>>>
>>>>library(zoo)
>>>>lst1<- split(d2,list(d$m1,d$n1)) 
>>>>d2<-do.call(rbind,lapply(lst1[lapply(lst1,nrow)!=0],function(x){ 
>>>>x[,9:14]<-NA; 
>>>>x[,9:10][x$Qm<=c11,]<-cumsum(x[,5:6][x$Qm<=c11,]); 
>>>>x[,11:12][x$Qn<=c12,]<-cumsum(x[,5:6][x$Qn<=c12,]);
>>>>x[,13:14]<-cumsum(x[,5:6]); 
>>>>colnames(x)[9:14]<- c("cterm1_P0L","cterm1_P1L","cterm1_P0H","cterm1_P1H","sumTerm1_p0","sumTerm1_p1"); 
>>>>x1<-na.locf(x); 
>>>>x1[,9:14][is.na(x1[,9:14])]<-0; 
>>>>x1}
>>>>)) 
>>>>row.names(d2)<-1:nrow(d2)
>>>>
>>>>head(d2)
>>>>tail(d2)
>>>>d2
>>>>
>>>>
>>>>res1<-aggregate(.~m1+n1,data=d2[,c(1:2,9:12)],max)
>>>>res1
>>>>
>>>>d3<-res1[res1[,4]<=beta & res1[,6]<beta,] 
>>>>tail(d3)
>>>>d3
>>>>
>>>>
>>>>if (nrow(d3)==0) break
>>>>
>>>>library(plyr) 
>>>>d4<- join(d3,d2,by=c("m1","n1"),type="inner")
>>>>head(d4) 
>>>>tail(d4)
>>>>#########################
>>>>################ 2nd stage ################################
>>>>#### this method not look at the combination of m1 and n1,
>>>>#### so need to merger the orginal data 'd2
>>>>
>>>>res2<-do.call(rbind,lapply(unique(d4$m1),function(m1) 
>>>>do.call(rbind,lapply(unique(d4$n1),function(n1)
>>>>do.call(rbind,lapply(unique(d4$x1),function(x1)
>>>>do.call(rbind,lapply(unique(d4$y1),function(y1)
>>>>
>>>>#do.call(rbind,lapply(0:m1,function(x1) 
>>>>#do.call(rbind,lapply(0:n1,function(y1) 
>>>>do.call(rbind,lapply((m1+2):(maxN-2-n1),function(m) 
>>>>do.call(rbind,lapply((n1+2):(maxN-m),function(n) 
>>>>do.call(rbind,lapply(x1:(x1+m-m1), function(x) 
>>>>do.call(rbind,lapply(y1:(y1+n-n1), function(y)
>>>>expand.grid(m1,n1,x1,y1,m,n,x,y)) ))))))))))))))) 
>>>>
>>>>names(res2)<- c("m1","n1","x1","y1","m","n","x","y") 
>>>>attr(res2,"out.attrs")<-NULL 
>>>>tail(res2)
>>>>
>>>>#install.packages(pkgs="plyr")
>>>>library(plyr) 
>>>>
>>>>res3<- join(res2,d4,by=c("m1","n1","x1","y1"),type="inner") 
>>>>head(res3)
>>>>tail(res3)
>>>>
>>>>res3<-res3[,c(1:16)]
>>>>head(res3)
>>>>tail(res3)
>>>>
>>>>############################################################# add more variables another method #############################
>>>>
>>>>res3<- within(res3,{
>>>>term2_p0<- dbinom(x1,m1, p0L, log=FALSE)* dbinom(y1,n1,p0H, log=FALSE)*dbinom(x-x1,m-m1, p0L, log=FALSE)* dbinom(y-y1,n-n1,p0H, log=FALSE);
>>>>term2_p1<- dbinom(x1,m1, p1L, log=FALSE)* dbinom(y1,n1,p1H, log=FALSE)*dbinom(x-x1,m-m1, p1L, log=FALSE)* dbinom(y-y1,n-n1,p1H, log=FALSE)}) 
>>>>
>>>>#term2_p0<- term1_p0*dbinom(x-x1,m-m1, p0L, log=FALSE)* dbinom(y-y1,n-n1,p0H, log=FALSE);
>>>>#term2_p1<- term1_p1*dbinom(x-x1,m-m1, p1L, log=FALSE)* dbinom(y-y1,n-n1,p1H, log=FALSE)}) 
>>>>
>>>>res4<-do.call(rbind,lapply(seq_len(nrow(res3)),function(i){
>>>>Pm2<-rbeta(1000,0.2+res3[i,"x"],0.8+res3[i,"m"]-res3[i,"x"]);
>>>>Pn2<- rbeta(1000,0.2+res3[i,"y"],0.8+res3[i,"n"]-res3[i,"y"]);
>>>>Fm2<- ecdf(Pm2); 
>>>>Fn2<- ecdf(Pn2);
>>>>#Fmm2<- Fm2(res3[i,"p1"]);
>>>>#Fnn2<- Fn2(res3[i,"p2"]);
>>>>
>>>>Fmm2<- Fm2(p1L);
>>>>Fnn2<- Fn2(p1H);
>>>>
>>>>R2<- (Fmm2+Fnn2)/2; 
>>>>Fmm_f2<- max(R2, Fmm2);
>>>>Fnn_f2<- min(R2, Fnn2); 
>>>>Qm2<- 1-Fmm_f2;
>>>>Qn2<- 1-Fnn_f2;
>>>>data.frame(Qm2,Qn2)}))
>>>>
>>>>res5<-cbind(res3,res4)
>>>>head(res5,5)
>>>>
>>>>################################################# calculate cumulative sum #################################
>>>>
>>>>lst2<- split(res5,list(res5$m1,res5$n1,res5$m,res5$n))
>>>>
>>>>res5<-do.call(rbind,lapply(lst2[lapply(lst2,nrow)!=0],
>>>>function(x){ 
>>>>x[,21:26]<-NA; 
>>>>x[,21:22][x$Qm2<=c1 & x$Qm>c11,]<-cumsum(x[,17:18][x$Qm2<=c1 & x$Qm>c11,]); 
>>>>x[,23:24][x$Qn2<=c2 & x$Qn>c12,]<-cumsum(x[,17:18][x$Qn2<=c2 & x$Qn>c12,]); 
>>>>x[,25:26]<-cumsum(x[,17:18]);
>>>>
>>>>colnames(x)[21:26]<- c("cterm2_P1L","cterm2_P0L","cterm2_P1H","cterm2_P0H","sumT2p0","sumT2p1"); 
>>>>x2<-na.locf(x); 
>>>>x2[,21:26][is.na(x2[,21:26])]<-0; x2}
>>>>)) 
>>>>row.names(res5)<-1:nrow(res5)
>>>>
>>>>head(res5)
>>>>tail(res5)
>>>>
>>>>res6<-aggregate(.~m1+n1+m+n,data=res5[,c(1:6,9:12,21:24)] ,max) 
>>>>
>>>>
>>>>res6<-within(res6,{AL<-1-(cterm1_P0L + cterm2_P0L);
>>>>AH<-1-(cterm1_P0H + cterm2_P0H);
>>>>BL<-(cterm1_P1L + cterm2_P1L);
>>>>BH<-(cterm1_P1H + cterm2_P1H);
>>>>EN<-m1 + (m-m1)*(1-cterm1_P0L) + n1 + (n-n1)*(1-cterm1_P0H);
>>>>N<-m+n
>>>>})
>>>>head(res6)
>>>>tail(res6)
>>>>
>>>>res7<-res6[,c(1:4,12,14,15:20)]
>>>>
>>>>res7<-res7[res7[,9]<=beta & res7[,10]<=beta &res7[,11] <= alpha & res7[,12] <= alpha,] 
>>>>res7
>>>>}
>>>>
>>>>
>>>>
>>>>
>>>>maxN<-9
>>>>p0L<-0.05
>>>>p0H<-0.05
>>>>p1L<-0.25
>>>>p1H<-0.25
>>>>
>>>>alpha<-0.4
>>>>beta<-0.4
>>>>
>>>>search(0.01,0.02,0.15,0.20)
>>>>
>>>>
>>>>
>>>>
>>>>On Sat, Mar 9, 2013 at 8:45 PM, arun <smartpink111 at yahoo.com> wrote:
>>>>
>>>>When I run the whole code:
>>>>>
>>>>>
>>>>> opt1
>>>>>[[1]]
>>>>>NULL
>>>>>
>>>>>[[2]]
>>>>>NULL
>>>>>
>>>>>[[3]]
>>>>>NULL
>>>>>
>>>>>[[4]]
>>>>>NULL
>>>>>
>>>>>[[5]]
>>>>>NULL
>>>>>
>>>>>[[6]]
>>>>>NULL
>>>>>
>>>>>[[7]]
>>>>>NULL
>>>>>
>>>>>[[8]]
>>>>>NULL
>>>>>
>>>>>[[9]]
>>>>>NULL
>>>>>
>>>>>[[10]]
>>>>>NULL
>>>>>
>>>>>> max1
>>>>>Error: object 'max1' not found
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>________________________________
>>>>>From: Joanna Zhang <zjoanna2013 at gmail.com>
>>>>>To: arun <smartpink111 at yahoo.com>
>>>>>Sent: Friday, March 8, 2013 3:07 PM
>>>>>
>>>>>Subject: Re: max row
>>>>>
>>>>>
>>>>>Arun,
>>>>>
>>>>>I still have some problems. I got a message that the computer was lack of memory as the data are too big. 
>>>>>So I wan to use run several loops to solve this problem. I need to get the 'opt1' and 'max1' for each combination of c11, c12, c1, c2, and loop through the different combinations of c11,c12,c1,c2, and output the opt1 and max1 along with the values of each combination. I include the code below, it's not working. Thanks very much in advance!
>>>>>
>>>>>maxN<-9
>>>>>
>>>>>p0L<-0.05
>>>>>p0H<-0.05
>>>>>p1L<-0.25
>>>>>p1H<-0.25
>>>>>
>>>>>alpha<-0.9
>>>>>beta<-0.9
>>>>>
>>>>>opt1<-vector('list', 10)
>>>>>max1<-vextor('list', 10)
>>>>>
>>>>>
>>>>>search<-function(c11,c12,c1,c2) {
>>>>>
>>>>>fopt<-matrix('list',(maxN-7),(maxN-6))
>>>>>fmax<-matrix('list',(maxN-7),(maxN-6))
>>>>>
>>>>>for (a in 2: (maxN-6)) {
>>>>>for (b in 2: (maxN-a-4)) {
>>>>>d <- data.frame ()
>>>>>for ( m1 in a:a) {
>>>>>     for (n1 in b: b){
>>>>>          for (x1 in 0: m1) {
>>>>>               for (y1 in 0: n1) {                                              
>>>>>
>>>>>                       term1_p0 = dbinom(x1,m1, p0L, log=FALSE)* dbinom(y1,n1,p0H, log=FALSE)
>>>>>                       term1_p1 = dbinom(x1,m1, p1L, log=FALSE)* dbinom(y1,n1,p1H, log=FALSE)                                                         
>>>>>
>>>>>                       d<-rbind(d, c(m1,n1,x1,y1,term1_p0,term1_p1))
>>>>>}}
>>>>>}}
>>>>>colnames(d)<-c("m1","n1","x1","y1","term1_p0","term1_p1") 
>>>>>#tail(d)
>>>>>
>>>>>########## add Qm Qn ##################################
>>>>>set.seed(8)
>>>>>d1<-do.call(rbind,lapply(seq_len(nrow(d)),function(i){
>>>>>Pm<- rbeta(1000,0.2+d[i,"x1"],0.8+d[i,"m1"]-d[i,"x1"]);
>>>>>Pn<- rbeta(1000,0.2+d[i,"y1"],0.8+d[i,"n1"]-d[i,"y1"]);
>>>>>Fm<- ecdf(Pm);
>>>>>Fn<- ecdf(Pn);
>>>>>
>>>>>Fmm<- Fm(p1L);
>>>>>Fnn<- Fn(p1H);
>>>>>
>>>>>R<- (Fmm+Fnn)/2; 
>>>>>Fmm_f<- max(R, Fmm);
>>>>>Fnn_f<- min(R, Fnn);
>>>>>Qm<- 1-Fmm_f;
>>>>>Qn<- 1-Fnn_f;
>>>>>data.frame(Qm,Qn)}))
>>>>>d2<-cbind(d,d1)
>>>>>#head(d2)
>>>>>
>>>>>
>>>>>library(zoo)
>>>>>lst1<- split(d2,list(d$m1,d$n1)) 
>>>>>d2<-do.call(rbind,lapply(lst1[lapply(lst1,nrow)!=0],function(x){ 
>>>>>x[,9:14]<-NA; 
>>>>>x[,9:10][x$Qm<=c11,]<-cumsum(x[,5:6][x$Qm<=c11,]); 
>>>>>x[,11:12][x$Qn<=c12,]<-cumsum(x[,5:6][x$Qn<=c12,]);
>>>>>x[,13:14]<-cumsum(x[,5:6]); 
>>>>>colnames(x)[9:14]<- c("cterm1_P0L","cterm1_P1L","cterm1_P0H","cterm1_P1H","sumTerm1_p0","sumTerm1_p1"); 
>>>>>x1<-na.locf(x); 
>>>>>x1[,9:14][is.na(x1[,9:14])]<-0; 
>>>>>x1}
>>>>>)) 
>>>>>row.names(d2)<-1:nrow(d2)
>>>>>
>>>>>#head(d2)
>>>>>#tail(d2)
>>>>>
>>>>>res1<-aggregate(.~m1+n1,data=d2[,c(1:2,9:12)],max)
>>>>>#head(res1)
>>>>>
>>>>>d3<-res1[res1[,4]<=beta & res1[,6]<beta,] 
>>>>>#tail(d3)
>>>>>
>>>>>if (nrow(d3)==0) break
>>>>>
>>>>>library(plyr) 
>>>>>d4<- join(d3,d2,by=c("m1","n1"),type="inner")
>>>>>#head(d4) 
>>>>>#tail(d4)
>>>>>
>>>>>
>>>>>
>>>>>res2<-do.call(rbind,lapply(unique(d4$m1),function(m1) 
>>>>>do.call(rbind,lapply(unique(d4$n1),function(n1)
>>>>>do.call(rbind,lapply(unique(d4$x1),function(x1)
>>>>>do.call(rbind,lapply(unique(d4$y1),function(y1)
>>>>>
>>>>>#do.call(rbind,lapply(0:m1,function(x1) 
>>>>>#do.call(rbind,lapply(0:n1,function(y1) 
>>>>>do.call(rbind,lapply((m1+2):(maxN-2-n1),function(m) 
>>>>>do.call(rbind,lapply((n1+2):(maxN-m),function(n) 
>>>>>do.call(rbind,lapply(x1:(x1+m-m1), function(x) 
>>>>>do.call(rbind,lapply(y1:(y1+n-n1), function(y)
>>>>>expand.grid(m1,n1,x1,y1,m,n,x,y)) ))))))))))))))) 
>>>>>
>>>>>names(res2)<- c("m1","n1","x1","y1","m","n","x","y") 
>>>>>attr(res2,"out.attrs")<-NULL 
>>>>>#tail(res2)
>>>>>
>>>>>#install.packages(pkgs="plyr")
>>>>>library(plyr) 
>>>>>
>>>>>res3<- join(res2,d4,by=c("m1","n1","x1","y1"),type="inner") 
>>>>>#head(res3)
>>>>>#tail(res3)
>>>>>
>>>>>res3<-res3[,c(1:16)]
>>>>>#head(res3)
>>>>>#tail(res3)
>>>>>
>>>>>############################################################# add more variables another method #############################
>>>>>
>>>>>res3<- within(res3,{
>>>>>term2_p0<- dbinom(x1,m1, p0L, log=FALSE)* dbinom(y1,n1,p0H, log=FALSE)*dbinom(x-x1,m-m1, p0L, log=FALSE)* dbinom(y-y1,n-n1,p0H, log=FALSE);
>>>>>term2_p1<- dbinom(x1,m1, p1L, log=FALSE)* dbinom(y1,n1,p1H, log=FALSE)*dbinom(x-x1,m-m1, p1L, log=FALSE)* dbinom(y-y1,n-n1,p1H, log=FALSE)}) 
>>>>>
>>>>>#term2_p0<- term1_p0*dbinom(x-x1,m-m1, p0L, log=FALSE)* dbinom(y-y1,n-n1,p0H, log=FALSE);
>>>>>#term2_p1<- term1_p1*dbinom(x-x1,m-m1, p1L, log=FALSE)* dbinom(y-y1,n-n1,p1H, log=FALSE)}) 
>>>>>
>>>>>
>>>>>res4<-do.call(rbind,lapply(seq_len(nrow(res3)),function(i){
>>>>>Pm2<-rbeta(1000,0.2+res3[i,"x"],0.8+res3[i,"m"]-res3[i,"x"]);
>>>>>Pn2<- rbeta(1000,0.2+res3[i,"y"],0.8+res3[i,"n"]-res3[i,"y"]);
>>>>>Fm2<- ecdf(Pm2); 
>>>>>Fn2<- ecdf(Pn2);
>>>>>#Fmm2<- Fm2(res3[i,"p1"]);
>>>>>#Fnn2<- Fn2(res3[i,"p2"]);
>>>>>
>>>>>Fmm2<- Fm2(p1L);
>>>>>Fnn2<- Fn2(p1H);
>>>>>
>>>>>R2<- (Fmm2+Fnn2)/2; 
>>>>>Fmm_f2<- max(R2, Fmm2);
>>>>>Fnn_f2<- min(R2, Fnn2); 
>>>>>Qm2<- 1-Fmm_f2;
>>>>>Qn2<- 1-Fnn_f2;
>>>>>data.frame(Qm2,Qn2)}))
>>>>>
>>>>>res5<-cbind(res3,res4)
>>>>>#head(res5,5)
>>>>>
>>>>>################################################# calculate cumulative sum #################################
>>>>>
>>>>>lst2<- split(res5,list(res5$m1,res5$n1,res5$m,res5$n))
>>>>>
>>>>>res5<-do.call(rbind,lapply(lst2[lapply(lst2,nrow)!=0],
>>>>>function(x){ 
>>>>>x[,21:26]<-NA; 
>>>>>x[,21:22][x$Qm2<=c1 & x$Qm>c11,]<-cumsum(x[,17:18][x$Qm2<=c1 & x$Qm>c11,]); 
>>>>>x[,23:24][x$Qn2<=c2 & x$Qn>c12,]<-cumsum(x[,17:18][x$Qn2<=c2 & x$Qn>c12,]); 
>>>>>x[,25:26]<-cumsum(x[,17:18]);
>>>>>
>>>>>colnames(x)[21:26]<- c("cterm2_P1L","cterm2_P0L","cterm2_P1H","cterm2_P0H","sumT2p0","sumT2p1"); 
>>>>>x2<-na.locf(x); 
>>>>>x2[,21:26][is.na(x2[,21:26])]<-0; x2}
>>>>>)) 
>>>>>row.names(res5)<-1:nrow(res5)
>>>>>
>>>>>#head(res5)
>>>>>#tail(res5)
>>>>>
>>>>>res6<-aggregate(.~m1+n1+m+n,data=res5,max) 
>>>>>#res6
>>>>>
>>>>>res6<-within(res6,{AL<-1-(cterm1_P0L + cterm2_P0L);
>>>>>AH<-1-(cterm1_P0H + cterm2_P0H);
>>>>>BL<-(cterm1_P1L + cterm2_P1L);
>>>>>BH<-(cterm1_P1H + cterm2_P1H);
>>>>>EN<-m1 + (m-m1)*(1-cterm1_P0L) + n1 + (n-n1)*(1-cterm1_P0H);
>>>>>N<-m+n
>>>>>})
>>>>>#head(res6)
>>>>>#tail(res6)
>>>>>
>>>>>res7<-aggregate(.~m1+n1+m+n,data=res6[,c(1:2,5:6,27:32)],max) 
>>>>>res7<-res6[,c(1:4,9,11,27:32)]
>>>>>
>>>>>res7<-res7[res7[,9]<=beta & res7[,10]<=beta &res7[,11] <= alpha & res7[,12] <= alpha,] 
>>>>>
>>>>>fopt[[a,b]]<-res7[which.min(res7$EN),] 
>>>>>fmax[[a,b]]<-res7[which.min(res7$N),]
>>>>>
>>>>>}}
>>>>>
>>>>>opt<-do.call(rbind,lapply(fopt,function(x) x[which.min(x$EN),]))
>>>>>max<-do.call(rbind,lapply(fmax,function(x) x[which.min(x$N),]))
>>>>>
>>>>>opt1<-opt[which.min(opt$EN),]
>>>>>max1<-max[which.min(max$N),]
>>>>>
>>>>>}
>>>>>opt1
>>>>>max1
>>>>>
>>>>>
>>>>>
>>>>>### loop c11,c12,c1, c2
>>>>>
>>>>>for (c11 in seq(0.05,0.06, by=0.01)) {
>>>>>for (c12 in seq(0.05,0.06, by=0.01)) {
>>>>>for (c1 in seq(0.2,0.3,by=0.1)){
>>>>>for (c2 in seq(0.2,0.3,by=0.1)){
>>>>>search(c11,c12,c1,c2) # create opt1 and max1
>>>>>par<-c(c11,c12,c1,c2);
>>>>>print(par)
>>>>>print(opt1)
>>>>>print(max1)
>>>>>
>>>>>}}}}
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>On Thu, Mar 7, 2013 at 10:54 AM, arun <smartpink111 at yahoo.com> wrote:
>>>>>
>>>>>No problem.
>>>>>>good luck.
>>>>>>Arun
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>________________________________
>>>>>> From: Joanna Zhang <zjoanna2013 at gmail.com>
>>>>>>To: arun <smartpink111 at yahoo.com>
>>>>>>Sent: Thursday, March 7, 2013 11:36 AM
>>>>>>Subject: Re: max row
>>>>>>
>>>>>>
>>>>>>
>>>>>>Thanks! I got it and solved the problem.
>>>>>>
>>>>>>
>>>>>>
>>>>>>On Wed, Mar 6, 2013 at 4:20 PM, arun <smartpink111 at yahoo.com> wrote:
>>>>>>
>>>>>>aggregate() gives the maximum values for each columns based on the combination.  Not the entire row.  Because, as you can see in a particular row, all the columns doesn't have the maximum values.
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>----- Original Message -----
>>>>>>>From: arun <smartpink111 at yahoo.com>
>>>>>>>To: Joanna Zhang <zjoanna2013 at gmail.com>
>>>>>>>Cc:
>>>>>>>Sent: Wednesday, March 6, 2013 5:11 PM
>>>>>>>Subject: Re: max row
>>>>>>>
>>>>>>>HI,
>>>>>>>
>>>>>>> res6<-aggregate(.~m1+n1+m+n,data=res5[,c(1:2,5:6,27:32)],max)
>>>>>>> res6
>>>>>>>#  m1 n1 m n N       EN        BH        BL        AH       AL
>>>>>>>#1  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1854938 0.097500
>>>>>>>#2  2  2 5 4 9 4.663488 0.3664627 0.7207031 0.1854938 0.097500
>>>>>>>#3  3  2 5 4 9 5.737688 0.3123894 0.6591797 0.2262191 0.142625
>>>>>>>#4  2  2 4 5 9 4.751481 0.4165192 0.6125565 0.1854938 0.097500
>>>>>>>#5  2  3 4 5 9 5.647438 0.3624458 0.6125565 0.2262191 0.097500
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>res7<-res5[,c(1:2,5:6,27:32)]
>>>>>>>res8<-split(res7, list(res7$m1,res7$n1,res7$m,res7$n))
>>>>>>>res9<-res8[lapply(res8,nrow)!=0]
>>>>>>>
>>>>>>> res9[[1]]
>>>>>>>   m1 n1 m n N       EN        BH        BL        AH         AL
>>>>>>>1   2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>2   2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>3   2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>4   2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>5   2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>6   2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>7   2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>8   2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>9   2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>10  2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>11  2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>12  2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>13  2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>14  2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>15  2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>16  2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>17  2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>18  2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>19  2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>20  2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>21  2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>22  2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>23  2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>24  2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>25  2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>26  2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>27  2  2 4 4 8 4.565988 0.3164062 0.5625000 0.1854938 0.09750000
>>>>>>>28  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>29  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>30  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>31  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>32  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>33  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>34  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>35  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>36  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>37  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>38  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>39  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>40  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>41  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>42  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>43  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>44  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>45  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>46  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>47  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>48  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>49  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>50  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>51  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>52  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>53  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>54  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>55  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>56  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>57  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>58  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>59  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>60  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>61  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>62  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>63  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>64  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>65  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>66  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>67  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>68  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>69  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>70  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>71  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>72  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>73  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>74  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>75  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>76  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>77  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>78  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>79  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>80  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>81  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1156600 0.02766627
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>Wha's wrong with row 81??
>>>>>>>
>>>>>>>Based on this, aren't you getting the maximum values for first combination???
>>>>>>>
>>>>>>>#  m1 n1 m n N       EN        BH        BL        AH       AL
>>>>>>>#1  2  2 4 4 8 4.565988 0.3831482 0.6292419 0.1854938 0.097500
>>>>>>>A.K.
>>>>>>>________________________________
>>>>>>>From: Joanna Zhang <zjoanna2013 at gmail.com>
>>>>>>>To: arun <smartpink111 at yahoo.com>
>>>>>>>Sent: Wednesday, March 6, 2013 4:01 PM
>>>>>>>Subject: max row
>>>>>>>
>>>>>>>
>>>>>>>Hi,
>>>>>>>
>>>>>>>I used aggregate to exact the max row of each block (combination  of m1, n1, m, n), but values of  AH, AL for m1 = 2, n1=2, m=4, n=4 (row 81) are not correct, they are the values of AH, AL for row 82. 
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>maxN<-9
>>>>>>>c11<-0.20
>>>>>>>c12<-0.20
>>>>>>>c1<-0.35
>>>>>>>c2<-0.35
>>>>>>>
>>>>>>>p0L<-0.05
>>>>>>>p0H<-0.05
>>>>>>>p1L<-0.25
>>>>>>>p1H<-0.25
>>>>>>>
>>>>>>>alpha<-0.60
>>>>>>>beta<-0.60
>>>>>>>
>>>>>>>d <- data.frame ()
>>>>>>>for ( m1 in 2: (maxN-6)) {
>>>>>>>     for (n1 in 2: (maxN-m1-4)){
>>>>>>>          for (x1 in 0: m1) {
>>>>>>>               for (y1 in 0: n1) {
>>>>>>>
>>>>>>>                       term1_p0 = dbinom(x1,m1, p0L, log=FALSE)* dbinom(y1,n1,p0H, log=FALSE)
>>>>>>>                       term1_p1 = dbinom(x1,m1, p1L, log=FALSE)* dbinom(y1,n1,p1H, log=FALSE)                                                         
>>>>>>>
>>>>>>>                       d<-rbind(d, c(m1,n1,x1,y1,term1_p0,term1_p1))
>>>>>>>}}
>>>>>>>}}
>>>>>>>colnames(d)<-c("m1","n1","x1","y1","term1_p0","term1_p1") 
>>>>>>>tail(d)
>>>>>>>
>>>>>>>set.seed(8)
>>>>>>>d1<-do.call(rbind,lapply(seq_len(nrow(d)),function(i){
>>>>>>>Pm<- rbeta(1000,0.2+d[i,"x1"],0.8+d[i,"m1"]-d[i,"x1"]);
>>>>>>>Pn<- rbeta(1000,0.2+d[i,"y1"],0.8+d[i,"n1"]-d[i,"y1"]);
>>>>>>>Fm<- ecdf(Pm);
>>>>>>>Fn<- ecdf(Pn);
>>>>>>>#Fmm<- Fm(d[i,"p11"]);
>>>>>>>#Fnn<- Fn(d[i,"p12"]);
>>>>>>>Fmm<- Fm(p1L);
>>>>>>>Fnn<- Fn(p1H);
>>>>>>>R<- (Fmm+Fnn)/2; 
>>>>>>>Fmm_f<- max(R, Fmm);
>>>>>>>Fnn_f<- min(R, Fnn);
>>>>>>>Qm<- 1-Fmm_f;
>>>>>>>Qn<- 1-Fnn_f;
>>>>>>>data.frame(Qm,Qn)}))
>>>>>>>d2<-cbind(d,d1)
>>>>>>>head(d2)
>>>>>>>
>>>>>>>
>>>>>>>library(zoo)
>>>>>>>lst1<- split(d2,list(d$m1,d$n1)) 
>>>>>>>d2<-do.call(rbind,lapply(lst1[lapply(lst1,nrow)!=0],function(x){ 
>>>>>>>x[,9:14]<-NA; 
>>>>>>>x[,9:10][x$Qm<=c11,]<-cumsum(x[,5:6][x$Qm<=c11,]); 
>>>>>>>x[,11:12][x$Qn<=c12,]<-cumsum(x[,5:6][x$Qn<=c12,]);
>>>>>>>x[,13:14]<-cumsum(x[,5:6]); 
>>>>>>>colnames(x)[9:14]<- c("cterm1_P0L","cterm1_P1L","cterm1_P0H","cterm1_P1H","sumTerm1_p0","sumTerm1_p1"); 
>>>>>>>x1<-na.locf(x); 
>>>>>>>x1[,9:14][is.na(x1[,9:14])]<-0; 
>>>>>>>x1}
>>>>>>>)) 
>>>>>>>row.names(d2)<-1:nrow(d2)
>>>>>>>
>>>>>>>head(d2)
>>>>>>>tail(d2)
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>res1<-aggregate(.~m1+n1,data=d2[,c(1:2,9:12)],max)
>>>>>>>
>>>>>>>d3<-res1[res1[,4]<=beta & res1[,6]<beta,] 
>>>>>>>tail(d3)
>>>>>>>d3
>>>>>>>
>>>>>>>
>>>>>>>if (nrow(d3)==0) break
>>>>>>>
>>>>>>>library(plyr) 
>>>>>>>d4<- join(d3,d2,by=c("m1","n1"),type="inner")
>>>>>>>head(d4) 
>>>>>>>tail(d4)
>>>>>>>
>>>>>>>res2<-do.call(rbind,lapply(unique(d4$m1),function(m1) 
>>>>>>>do.call(rbind,lapply(unique(d4$n1),function(n1)
>>>>>>>do.call(rbind,lapply(unique(d4$x1),function(x1)
>>>>>>>do.call(rbind,lapply(unique(d4$y1),function(y1)
>>>>>>>
>>>>>>>do.call(rbind,lapply((m1+2):(maxN-2-n1),function(m) 
>>>>>>>do.call(rbind,lapply((n1+2):(maxN-m),function(n) 
>>>>>>>do.call(rbind,lapply(x1:(x1+m-m1), function(x) 
>>>>>>>do.call(rbind,lapply(y1:(y1+n-n1), function(y)
>>>>>>>expand.grid(m1,n1,x1,y1,m,n,x,y)) ))))))))))))))) 
>>>>>>>
>>>>>>>names(res2)<- c("m1","n1","x1","y1","m","n","x","y") 
>>>>>>>attr(res2,"out.attrs")<-NULL 
>>>>>>>tail(res2)
>>>>>>>
>>>>>>>#install.packages(pkgs="plyr")
>>>>>>>library(plyr) 
>>>>>>>
>>>>>>>res3<- join(res2,d4,by=c("m1","n1","x1","y1"),type="inner") 
>>>>>>>head(res3)
>>>>>>>tail(res3)
>>>>>>>
>>>>>>>res3<-res3[,c(1:16)]
>>>>>>>head(res3)
>>>>>>>tail(res3)
>>>>>>>
>>>>>>>############################################################# add more variables another method #############################
>>>>>>>
>>>>>>>res3<- within(res3,{
>>>>>>>term2_p0<- dbinom(x1,m1, p0L, log=FALSE)* dbinom(y1,n1,p0H, log=FALSE)*dbinom(x-x1,m-m1, p0L, log=FALSE)* dbinom(y-y1,n-n1,p0H, log=FALSE);
>>>>>>>term2_p1<- dbinom(x1,m1, p1L, log=FALSE)* dbinom(y1,n1,p1H, log=FALSE)*dbinom(x-x1,m-m1, p1L, log=FALSE)* dbinom(y-y1,n-n1,p1H, log=FALSE)}) 
>>>>>>>
>>>>>>>
>>>>>>>res4<-do.call(rbind,lapply(seq_len(nrow(res3)),function(i){
>>>>>>>Pm2<-rbeta(1000,0.2+res3[i,"x"],0.8+res3[i,"m"]-res3[i,"x"]);
>>>>>>>Pn2<- rbeta(1000,0.2+res3[i,"y"],0.8+res3[i,"n"]-res3[i,"y"]);
>>>>>>>Fm2<- ecdf(Pm2); 
>>>>>>>Fn2<- ecdf(Pn2);
>>>>>>>
>>>>>>>Fmm2<- Fm2(p1L);
>>>>>>>Fnn2<- Fn2(p1H);
>>>>>>>
>>>>>>>R2<- (Fmm2+Fnn2)/2; 
>>>>>>>Fmm_f2<- max(R2, Fmm2);
>>>>>>>Fnn_f2<- min(R2, Fnn2); 
>>>>>>>Qm2<- 1-Fmm_f2;
>>>>>>>Qn2<- 1-Fnn_f2;
>>>>>>>data.frame(Qm2,Qn2)}))
>>>>>>>
>>>>>>>res5<-cbind(res3,res4)
>>>>>>>head(res5,5)
>>>>>>>
>>>>>>>################################################# calculate cumulative sum #################################
>>>>>>>
>>>>>>>lst2<- split(res5,list(res5$m1,res5$n1,res5$m,res5$n))
>>>>>>>
>>>>>>>res5<-do.call(rbind,lapply(lst2[lapply(lst2,nrow)!=0],
>>>>>>>function(x){ 
>>>>>>>x[,21:26]<-NA; 
>>>>>>>x[,21:22][x$Qm2<=c1 & x$Qm>c11,]<-cumsum(x[,17:18][x$Qm2<=c1 & x$Qm>c11,]); 
>>>>>>>x[,23:24][x$Qn2<=c2 & x$Qn>c12,]<-cumsum(x[,17:18][x$Qn2<=c2 & x$Qn>c12,]); 
>>>>>>>x[,25:26]<-cumsum(x[,17:18]);
>>>>>>>
>>>>>>>colnames(x)[21:26]<- c("cterm2_P1L","cterm2_P0L","cterm2_P1H","cterm2_P0H","sumT2p0","sumT2p1"); 
>>>>>>>x2<-na.locf(x); 
>>>>>>>x2[,21:26][is.na(x2[,21:26])]<-0; x2}
>>>>>>>)) 
>>>>>>>row.names(res5)<-1:nrow(res5)
>>>>>>>
>>>>>>>head(res5)
>>>>>>>tail(res5)
>>>>>>>
>>>>>>>
>>>>>>>res5<-within(res5,{AL<-1-(cterm1_P0L + cterm2_P0L);
>>>>>>>AH<-1-(cterm1_P0H + cterm2_P0H);
>>>>>>>BL<-(cterm1_P1L + cterm2_P1L);
>>>>>>>BH<-(cterm1_P1H + cterm2_P1H);
>>>>>>>EN<-m1 + (m-m1)*(1-cterm1_P0L) + n1 + (n-n1)*(1-cterm1_P0H);
>>>>>>>N<-m+n
>>>>>>>})
>>>>>>>head(res5)
>>>>>>>tail(res5)
>>>>>>>
>>>>>>>res6<-aggregate(.~m1+n1+m+n,data=res5[,c(1:2,5:6,27:32)],max) 
>>>>>>>res6                                                   
>>>>>>>
>>>>>>>
>>>>>>
>>>>>>
>>>>>>                                                                  
>>>>>  
>>>>
>>>
>>
>



More information about the R-help mailing list