R code for Formulas 1, 2, 3, and 4 in the paper
Recursive Formulas for Multinomial Probabilities with Applications


FORMULA 1

n<-150
k<-8
p<-c(0.1,0.1,0.1,0.1,0.1,0.1,0.2,0.2)
l<-c(8,8,8,8,8,8,8,8)
u<-c(30,30,30,30,30,30,30,30)

l1<-rev(cumsum(rev(l)))
l2<-cumsum(l)
u1<-rev(cumsum(rev(u)))
u2<-cumsum(u)

g<-seq(0,n+1,1)
g1<-seq(0,n+1,1)
for (i in 0:n){g[i+1]<-1/factorial(n-i)}

for (i in (k-2):1){
for (w in 0:n){b<-0
for (j in max(n-u1[i+2],l[i+1]+w):min(n-l1[i+2],u[i+1]+w)){
b<-b+(p[i+1]/p[i+2])^j*g[j+1]/factorial(max(0,j-w))}
g1[w+1]<-b}
g<-g1
g1<-seq(0,n+1,1)}

A<-0
for (i3 in l[1]:u[1]){A<-A+(p[1]/p[2])^i3*g[i3+1]/factorial(i3)}
A<-A*factorial(n)*(p[k]^n)
A





FORMULA 2

n<-100
k<-8
p<-c(1,1,1,1,1,1,1,2)/9

p1 <- p[-k]/(1-p[k])  
n1<- ceiling((n-1)/k)+1
n2<- floor(n/2)+1

M1<-0
for (ii in n1:(n2-1)){m<-max(0,n-(k-1)*ii+k-2)
kk<-k-1
nn<-n-ii
l<-rep(m,kk)
u<-rep(ii-1,kk)
l1<-rev(cumsum(rev(l)))
l2<-cumsum(l)
u1<-rev(cumsum(rev(u)))
u2<-cumsum(u)

g<-seq(0,nn+1,1)
g1<-seq(0,nn+1,1)
for (i in 0:nn){g[i+1]<-1/factorial(nn-i)}

for (i in (kk-2):1){
for (w in 0:nn){b<-0
for (j in max(nn-u1[i+2],l[i+1]+w):min(nn-l1[i+2],u[i+1]+w)){
b<-b+(p1[i+1]/p1[i+2])^j*g[j+1]/factorial(max(0,j-w))}
g1[w+1]<-b}
g<-g1
g1<-seq(0,nn+1,1)}

A<-0
for (i3 in l[1]:u[1]){A<-A+(p1[1]/p1[2])^i3*g[i3+1]/factorial(i3)}
A<-A*factorial(nn)*(p1[kk]^nn)
M1<-M1+A*dbinom(ii,n,p[k])}
M1<-M1+1-pbinom(n2-1,n,p[k])
M1




FORMULA 3

n<-100
k<-8
p<-c(1,1,1,1,1,1,1,2)/9

p1 <- p[-k]/(1-p[k])  
n1<- ceiling(n/k)
n2<- floor((n-1)/2)+1

M2<-0
for (ii in n1:(n2-1)){m<-max(0,n-(k-1)*ii)
kk<-k-1
nn<-n-ii
l<-rep(m,kk)
u<-rep(ii,kk)
l1<-rev(cumsum(rev(l)))
l2<-cumsum(l)
u1<-rev(cumsum(rev(u)))
u2<-cumsum(u)

g<-seq(0,nn+1,1)
g1<-seq(0,nn+1,1)
for (i in 0:nn){g[i+1]<-1/factorial(nn-i)}

for (i in (kk-2):1){
for (w in 0:nn){b<-0
for (j in max(nn-u1[i+2],l[i+1]+w):min(nn-l1[i+2],u[i+1]+w)){
b<-b+(p1[i+1]/p1[i+2])^j*g[j+1]/factorial(max(0,j-w))}
g1[w+1]<-b}
g<-g1
g1<-seq(0,nn+1,1)}

A<-0
for (i3 in l[1]:u[1]){A<-A+(p1[1]/p1[2])^i3*g[i3+1]/factorial(i3)}
A<-A*factorial(nn)*(p1[kk]^nn)
M2<-M2+A*dbinom(ii,n,p[k])}
M2<-M2+1-pbinom(n2-1,n,p[k])
M2




FORMULA 4 (k=3)

n<-80
k<-3
p<-c(1,2,3)/6

s1<-0
for (y1 in 0:floor(n/k)){s2<-0
for (y2 in (2*y1):(floor((n+(k-2)*y1)/(k-1)))){
s2<-s2+(p[2]/p[3])^y2/(factorial(max(0,y2-y1))*factorial(n-y2))}
s1<-s1+s2*(p[1]/p[2])^y1/factorial(y1)}
O<-s1*factorial(n)*p[k]^n
O






FORMULA 4 (k=4)

n<-80
k<-4
p<-c(1,2,3,4)/10

g<-matrix(0,n+2,n+2)
for (w in 0:((k-2)*n/k)){
g[w+2,w+2]<-(p[k-1]/p[k])^w/factorial(n-w)
for (z in (w+1):((n+w)/2)){
g[w+2,z+2]<-g[w+2,z+1]+(p[k-1]/p[k])^z/(factorial(z-w)*factorial(n-z))}}

s1<-0
for (y1 in 0:floor(n/k)){s2<-0
for (y2 in (2*y1):(floor((n+(k-2)*y1)/(k-1)))){
s2<-s2+(p[2]/p[3])^y2*(g[y2+2,floor((n+(k-3)*y2)/(k-2))+2]-g[y2+2,2*y2-y1+1])/factorial(max(0,y2-y1))}
s1<-s1+s2*(p[1]/p[2])^y1/factorial(y1)}
O<-s1*factorial(n)*p[k]^n
O





FORMULA 4 (k>4)

n<-160
k<-8
p<-c(1,2,3,4,5,6,7,8)/36

g<-matrix(0,n+2,n+2)
for (w in 0:((k-2)*n/k)){g[w+2,w+2]<-(p[k-1]/p[k])^w/factorial(n-w)
for (z in (w+1):((n+w)/2)){
g[w+2,z+2]<-g[w+2,z+1]+(p[k-1]/p[k])^z/(factorial(z-w)*factorial(n-z))}}

for (i in (k-4):1){g1<-matrix(0,n+2,n+2)
for (w in 0:((i+1)*n/k)){
for (z in (w):((n+(k-2-i)*w)/(k-1-i))){
g1[w+2,z+2]<-g1[w+2,z+1]+(p[i+2]/p[i+3])^z*(g[z+2,floor((n+(k-3-i)*z)/(k-2-i))+2]-g[z+2,2*z-w+1])/factorial(z-w)}}
g<-g1}

s1<-0
for (y1 in 0:floor(n/k)){s2<-0
for (y2 in (2*y1):(floor((n+(k-2)*y1)/(k-1)))){
s2<-s2+(p[2]/p[3])^y2*(g[y2+2,floor((n+(k-3)*y2)/(k-2))+2]-g[y2+2,2*y2-y1+1])/factorial(max(0,y2-y1))}
s1<-s1+s2*(p[1]/p[2])^y1/factorial(y1)}
O<-s1*factorial(n)*p[k]^n
O