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 <zjoanna2...@gmail.com>
To: arun <smartpink...@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 <smartpink...@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 <zjoanna2...@gmail.com>
>To: arun <smartpink...@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 <smartpink...@yahoo.com> wrote:
>
>I hope your problem is solved.
>>
>>
>>
>>
>>
>>
>>
>>________________________________
>>From: Joanna Zhang <zjoanna2...@gmail.com>
>>To: arun <smartpink...@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 <smartpink...@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 <zjoanna2...@gmail.com>
>>>To: arun <smartpink...@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 <smartpink...@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 <zjoanna2...@gmail.com>
>>>>To: arun <smartpink...@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 <smartpink...@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 <zjoanna2...@gmail.com>
>>>>>To: arun <smartpink...@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 <smartpink...@yahoo.com> wrote:
>>>>>
>>>>>No problem.
>>>>>>good luck.
>>>>>>Arun
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>________________________________
>>>>>> From: Joanna Zhang <zjoanna2...@gmail.com>
>>>>>>To: arun <smartpink...@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 <smartpink...@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 <smartpink...@yahoo.com>
>>>>>>>To: Joanna Zhang <zjoanna2...@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 <zjoanna2...@gmail.com>
>>>>>>>To: arun <smartpink...@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
>>>>>>>
>>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>
>>>
>>
>
______________________________________________
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.
