#Source for stratum program
#inputs:
# k:vector of stratum costs (length I)
# n:vector of stratum sizes (length I)
# s:vector of stratum variances (length I)
#
k<-c(2.5,4.2,6.7)
n<-c(3,5,7)
s<-c(7,9,11)
#outputs:
# ans: dynamically optimal permutation (length
# length(n))
# cost: vector of cumulative costs (length
# length(n))
# var: vector of variances (length length(n))
#the first job is to compute the start and end of each stratum . The end is easy. The start starts with 1, adds one to each cumulative sum, and lops off the last element of the vector.
start<-c(1,(cumsum(n)+1))[-(length(n)+1)]
end<-cumsum(n)
#the second job is to prepare x, a vector of length
# the number of items, ie length(n), which has
# a random permutation of each stratum, catenated.
x<-rep(0,sum(n))
for (i in 1:length(n)) x[start[i]:end[i]]_sample(start[i]:end[i])
#third, the vector del is prepared, which is the incremental decrease in variance. Note that the first element of each stratum requires special handling, since without it, the total variance would be infinite
b<-rep(s*(n**2),n)/(sum(n)**2)
m<-max(n)
w<-c(0,1/((1:m)*((1:m)+1)))
d<-rep(0,sum(n))
for (i in 1:length(n)) d[start[i]:end[i]]<-w[1:n[i]]
del<-b*d
# now del is prepared, with 0's at each "start" position
a<-del/rep(k,n)
a[start]<-(1+max(a))
#the last command ensures that each stratum is allocated one sample point before the del allocation begins
ord<-sort.list(-a)
#ord now has the optimal permutation
ans<-x[ord]
ans
cost<-cumsum(rep(k,n)[ord])[-(1:(length(n)-1))]
v<-rep(0,sum(n))
v[-start]<-(-del[-start])
v[start]<-((n**2)*s*(1-1/n))/((sum(n))**2)
var<-cumsum(v[ord])[-(1:(length(n)-1))]
#without lopping off the first I-1 elements of var, the variance would appear to climb before it starts to decrease. Actually it is infinite at those points. So I have taken those elemnts away, both on cost and on var.