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? ? ? ? ? ? ? ? ? ? ? ? ??? ? ? ? ? ? ? ? ? ? ? ???
>>>>>>>
>>>>>>>
>>>>>>
>>>>>>
>>>>>>? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ???
>>>>>??
>>>>
>>>
>>
>