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