> L<-function(rho,eta,xi,z,C,n,pn) + { + L<-0 + for(j in 1:nrow(n)) + { + Lfac<-1 + A<-length(z) + for(a in 1:A) + { + if(a==A) shape<-rho*(1-xi)*n[j,a] else + shape<-rho*((1-xi)*n[j,a]+xi*n[j,a+1]) + if(z[a] > fn<-function(par) + { + -log(L(par[1],par[2],par[3],z,C,n,pn)) + } > > C<-0.001 > > # alleles 20 to 24 inclusive, z are peak heights, alpha are Dirichlet-Multinomial parameters, C is cut-off > z<-c(0,100,0,400,0) # peak heights > alpha<-100*c(0.12748, 0.18543, 0.21854, 0.13411, 0.13576) # database numbers to match USCaucasian2, M=100 > > # possible \mathbf{n} values and their probabilities under DM model, hypothesis of 1 unknown actor > A<-length(z) > nn<-A*(A+1)/2 > n<-matrix(0,nn,A) > pn<-rep(0,nn) > k<-0 > for(i in 1:A) for(j in i:A) + { + k<-k+1 + n[k,i]<-n[k,i]+1 + n[k,j]<-n[k,j]+1 + pn[k]<-(1+(i!=j))*alpha[i]*alpha[j]+(i==j)*alpha[i] + } > pn<-pn/sum(pn) > > fit<-nlminb(c(2,100,0.1),fn,lower=c(0,0,0)) > fit$lik<-(-fit$obj/log(10)) > fit $par [1] 2.394166 104.420468 0.000000 $objective [1] 15.26963 $convergence [1] 0 $iterations [1] 19 $evaluations function gradient 20 74 $message [1] "relative convergence (4)" $lik [1] -6.631515 > > # Repeat for vanilla (multinomial) case > k<-0 > for(i in 1:A) for(j in i:A) + { + k<-k+1 + pn[k]<-(1+(i!=j))*alpha[i]*alpha[j] + } > pn<-pn/sum(pn) > > fit<-nlminb(c(2,100,0.1),fn,lower=c(0,0,0)) > fit$lik<-(-fit$obj/log(10)) > fit $par [1] 2.394166 104.420468 0.000000 $objective [1] 15.25722 $convergence [1] 0 $iterations [1] 19 $evaluations function gradient 20 74 $message [1] "relative convergence (4)" $lik [1] -6.626129