Is there an optimal way to find the best division of an interval of some positive integers? I am struggling with a conceptual problem.
I have positive integers from an interval [1800, 1850].
For every integer from that interval, let's say (without loss of generality) 1820, I have about 3000 horses. The 1820 number is a year of birth for a horse. Thoses horses were fed with a traditional food and some of those horses were fed with experimental food (there were 29 types of different experimental food). For every horse there was recorded a variable for each feeding named goodness of sneeze (the higer the goodness variable is, the better). Let's assume after every feeding a horse did sneeze. Every single horse could be fed with different type of food every time he came on feeding (with uniform distribution). Let us assume that sneeze for horses comes from Poisson distribution with lamba=1 parameter.
Now I am looking for the best [1800,1850] interval division on intervals like:
[1800,1810), [1810,1826), [1826,1850]

to say: for every subinterval this or that experimental food (or maybe traditional in some cases) gave best average sneeze for horses born in that interval. 
I do not know if it is needed, but let's assume that horses does not come on feeding with regularity. Some of them come more often than others. Experiment took 20 days.
If there is a good way of generating the best interval in a relatively fast way?
I tried to make a loop for i in 1 to 50 where i is a number of [1800,1850] interval divisions centers.
If i=1
I check:
[1800,1801],(1802,1850] 
[1800,1802],(1803,1850] 
...
[1800,1849],(1849,1850] 

and check which experimental food gave the biggest mean sneeze in that subinterval and answer the problem as this example:
[1800,1807], (1807,1850]

is the best division from division with 1 interval centers for horses born in [1800,1807] the best food is experimentalFoodnr25 and for horses born in (1807,1850] the best food is experimentalFoodnr14. 
With respect to traditional food they give 0,04 higher mean sneeze for horses. (0.04 is of course a weighted mean with respect to number of horses in both intervals)
Then I can go for i=2, and so on and so on but there higher the i is, the less horses are in the subintervals and the estimate of the average sneeze has greater standard error.
So I thought about to choose the best [1800,1850] division that has the biggest
weighted mean of a's where a is calculated from subinterval and is to be as formula:
$a = \phi( 1- p )^{-1} \times \sqrt{ Var(X)/n_{x} + Var(Y)/n_{y} } + \mu_{X} - \mu_{Y}$
where $X$ are the records for horses treated with the experimental food giving the highest average sneeze in that subinterval, $Y$ are the records for horses treated with traditional food in that subinterval. $\mu$ are means of that records, $Var$ are variances and p is the probability of that $P( \mu_{X}-\mu_{Y}>a)=p$  (where I assume $\mu_{X}$ has normal distributions) and $\phi$ is a standard normal distribution function and n's are number of records.
Can someone has any idea of relatively fast algorithm for that problem?
If the problem is not clear please tell me what to specify.
 A: I think that the most simple solution would include two steps:


*

*Find the best solution for even-sized clusters,

*Re-assign individual groups cluster by cluster and check if it results a better solution.


Below I post a code in R to give an example. I changed labels of your groups from 1800:1850 to 1:50 for simplicity and assumed that each group has a fixed size $n$. data is a vector sorted by groups (first 1:n cases is group 1, (n+1):(2*n+1) is a second group etc.). I used recursion since the number and labels of clusters can change (so for-loop was not an option).
assessClusters <- function(data, cluster, n)
                    tapply(data, rep(cluster, each=n), mean)

# re-assign a single value of given cluster

changeOne <- function(data, cluster, i, n) {
  crit <- assessClusters(data, cluster, n=n)

  # find the last value from the cluster

  k <- max(seq_along(cluster)[cluster == i])

  # make the cluster bigger

  if (k+1 <= length(cluster)) {
    cluster[k+1] <- cluster[k+1]-1
    if (mean(crit) < mean(assessClusters(data, cluster, n=n))) {
      cluster <- changeOne(data, cluster, i, n=n)
    } else {
      cluster[k+1] <- cluster[k+1]+1
    }
  }

  # make it smaller

  cluster[k] <- cluster[k]+1
  if (mean(crit) < mean(assessClusters(data, cluster, n=n))) {
    cluster <- changeOne(data, cluster, i, n=n)
  } else {
    cluster[k] <- cluster[k]-1
    return(cluster)
  } 
}

# iterate through all the clusters 

clustersIterate <- function(i=1, data, cluster, n) {
  cluster <- changeOne(data, cluster, i, n=n)
  if (i+1 < max(cluster)) {
    return(clustersIterate(i+1, data, cluster, n=n))
  } else {
    return(cluster)
  }
}

And here you can check how it works:
n <- 100
groups <- rep(1:50, each=n)
data <- rpois(50*n, 1)

# find the optimal number of clusters for even sized clusters
# max.even.clusters is the maximal number of clusters to check

max.even.clusters <- 10
tmp <- 0
for (i in 2:max.even.clusters) {
  tmp[i] <- mean(assessClusters(data,
                                rep(1:50, each=i, length.out=50),
                                n=n))
}
cluster <- rep(1:50, each=which.max(tmp), length.out=50)

# cls are the new labels for clusters

cls <- clustersIterate(1, data, cluster, n=n)

# below you can check if it actually works

assessClusters(data, cluster, n=n)
assessClusters(data, cls, n=n)

mean(assessClusters(data, cluster, n=n))
mean(assessClusters(data, cls, n=n))

Of course you can use different decision criteria then the mean (e.g. weighted mean, as in your question). Using mean could be tricky since it is not outlier-proof, however since the algorithm starts with even clusters and then re-assigns one group at a time, there should be no problems with outliers. The algorithm is certainly faster then looping through every possible solution. I don't argue that this is a best solution, however it points some possible directions. The code possibly contains bugs so if you'd like to use it then I would suggest making some tests.

Generally you should think of it as a clustering problem. The code above was inspired partly by k-means algorithm, however you can think also of a hierarchical clustering solution for it. The most common problem with cluster analysis is deciding on the number of clusters, but in your case it is simpler since you have criteria you want to maximize (cluster means). No matter what solution you'll choose you should take into consideration that there could be multiple cluster solutions for this problem, so you could end up in a local maxima. If you use the algorithm posted above, then what you can do is to start it with a different number of even-sized clusters and compare the solutions to check if the results are stable.
