> 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]<C) L1<-pgamma(C,shape,scale=eta) else L1<-dgamma(z[a],shape,scale=eta)
+ Lfac<-Lfac*L1
+ }
+ L<-L+pn[j]*Lfac
+ }
+ L
+ }
> 
> 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