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:

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.

  • $\begingroup$ Why do you want to find this optimal division? What phenomenon exactly are you trying to understand? I suspect that finding the optimal split as you described is going to be tough, but there might be other, easier ways to achieve the same end goal. $\endgroup$ – Ben Kuhn Feb 3 '15 at 21:22
  • $\begingroup$ Somehow this division will be an answer for preferences analysis. Isn't it? $\endgroup$ – Marcin Kosiński Feb 4 '15 at 21:57
  • $\begingroup$ I have no idea what you mean. When I look up "preferences analysis" I get this wiki page which makes no mention of time series of preferences at all. Can you elaborate? $\endgroup$ – Ben Kuhn Feb 4 '15 at 22:16
  • $\begingroup$ This does not mean the topic does not exist if it is not on a wiki. This is the problem to be solved. Does it matter from which field of science does it come? $\endgroup$ – Marcin Kosiński Feb 5 '15 at 10:24
  • $\begingroup$ I'm not claiming that the topic doesn't exist! Merely that it's not common knowledge what you mean by it and so the phrase "this will be an answer for preference analysis" doesn't mean anything to me. Anyway, hopefully you can find someone on this site who already knows how this is "an answer for preference analysis." Good luck! $\endgroup$ – Ben Kuhn Feb 5 '15 at 19:09

I think that the most simple solution would include two steps:

  1. Find the best solution for even-sized clusters,
  2. 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

# 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 {

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),
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.

  • $\begingroup$ Thanks @Tim ! I'll have a look at this in the evening. The idea of clustering is amazing. Great job! You've saved the day :) $\endgroup$ – Marcin Kosiński Feb 10 '15 at 17:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.