Approximating a Complementary Cumulative Distribution Function via a piece-wise function I hope this is not too much to read, but I tried to give you a specific overview over my problem.
I am currently trying to model the German electricity market, with a special focus on balancing power.
Considering this setting I have a data set consisting of hourly observations of balancing power in Megawatts for the year of 2013. As the application of this data set would not represent any uncertainty in the data I wanted to use a little trick to "assume" uncertainty. While this is really only a basic approximation my model has to have some sort of "uncertainty", else there will be a big mistake. In order to do this, I divided my data set into five blocks (too many more would be to much for my calculation) of varying Megawatt-steps calculating the "probability" of a specific value to be in one of these blocks using the Complementary cumulative distribution function.
What I mean by that is that I took the aforementioned function and tried to best approximate it in a (5-)piece-wise fashion. Unfortunately I had to do this in my example by a "sense of proportion". I hope the picture makes it a little more clear (On the x-axis the values are in Megawatts). 
EDIT: To be clear, I am not trying to best-possibly approximate the whole-function in a sense of a non parametric boundary estimation as user603 kindly described in his answer. (Or maybe this can be done using non parametric boundary estimation...?) I just want to fit the five (or whatever number) blocks best-possible to the function in a fashion that, e.g. the Sum of squared Errors is at a minimum. Also: i would love the knot location to be estimated.

The following code example should further explain what I do:
# This function divides the input data into 5 different blocks.
# v, w, x and y are all Megawatt values dividing the data set
getCall_block <- function(v, w, x, y, data){

  n <- nrow(data)
  #' parameters for function
  data_v <- subset(data, data[,2] <= v, select = c(colnames(data)[1],
                                                   colnames(data)[2]))
  data_w <- subset(data, data[,2] > v & data[,2] <= w, select = c(colnames(data)[1],
                                                                  colnames(data)[2]))
  data_x <- subset(data, data[,2] > w & data[,2] <= x, select = c(colnames(data)[1],
                                                                  colnames(data)[2]))
  data_y <- subset(data, data[,2] > x & data[,2] <= y, select = c(colnames(data)[1],
                                                                  colnames(data)[2]))
  data_r <- subset(data, data[,2] > y, select = c(colnames(data)[1],
                                                  colnames(data)[2]))
  #' 
  df_v <- getAverage(data, data_v, 0, v)
  df_w <- getAverage(data, data_w, v, w)
  df_x <- getAverage(data, data_x, w, x)
  df_y <- getAverage(data, data_y, x, y)
  df_r <- getAverage(data, data_r, y, max(data[,2]))

  return(rbind(df_v, df_w, df_x, df_y, df_r))
}

And the getAverage function...
#This function is called in the above mentioned function. 
#It calculates the mean average of the points of the CCDF lying in one block. 
#v and w are the input parameters in Megawatts. 
#They can be understood as starting point and end point of a given block.

getAverage <- function(data_c, data, v, w){
  #' Imposing greater equal on v if v > 0
  if(w > v){
    if (v > 0){
      v <- v + 1
    } else {
      v <- 0
    }
  #' Initializing div, df, mean, n, sum
    div <- 0
    mean <- 0
    n  <- nrow(data)
    sum <- 0
    df <- data.frame(ID = data[,1], call_freq = numeric(n))
  #' Creating inverse ECDF based on dataset
    F1 <- ecdf(data_c[,2])
  #' Building sum over every data point in the interval [v;w]/(v;w]
    for (i in seq(v, w, 1)){
      sum <- sum + (1 - F1(i))
      div <- div + 1
    }
    mean <- sum / div
    df$call_freq = mapply(function(x) return (x),  mean)
    return (merge(data, df, by = "ID"))
    } else {
      cat("Check input parameters. v is greater than w.")
  }
}

My questions are the following:
I) Does this way of approximating a CCDF sound okay for you guys? What did I overlook or completely mess up?
II) Is there a way to automatically approximate the block sizes, in order to not rely on a sense of proportion. E.g. minimizing the sum of squared errors between the piece-wise function and the original function. Maybe this could be done by using a clustering algorithm like k-means? But then how to cluster this CCDF...?
III) Is there a far more easier way to approximate any function in a piece-wise fashion?
Thanks in advance!
 A: First, It's not really clear to me what you are trying to do. Your question is indeed too long winded and I think this is preventing other users on our site from offering you more answers to choose from. 
Anyway, I will take a stab at it, hopefully it'll help restructure you question. In the convention, tools and approaches I will describe I will assume that the function you try to approximate is monotone increasing. In your case this means working with $1-\hat{F}$ (1 minus the empirical cdf) but that shouldn't be a problem.
Now, to answer your questions:

I) Does this way of approximating a CCDF sound okay for you guys? What
  did I overlook or completely mess up?

No: you should add more constraints to your optimization problem. For example, you know --through elementary subject knowledge-- that the CDF (CCDF) is monotonous increasing (decreasing). You should explicitely constraint the approximations to also be. 
Including this constraint will improve the finite sample performance of the estimates and will lessen the dependence of the result on the choice of the nkot points. As to how to do so, I explain in the answer to your question iii) below. Also, you probably want a lower (upper) bound approximation to the (CCDF) CDF for all $x<x_M$ where $x_M$ is your estimate of the largest possible observation. You should also add this as a constraint in the fitting procedure. 

II) Is there a way to automatically approximate the block sizes, in
  order to not rely on a sense of proportion. E.g. minimizing the sum of
  squared errors between the piece-wise function and the original
  function. Maybe this could be done by using a clustering algorithm
  like k-means? But then how to cluster this CCDF...?

If you consider the knot location as fixed, yes: this is a convex optimization problem. If you consider the knot locations as a parameter, then the resulting problem is not convex but can be approximated by a convex optimization problem. In principle, this could be done by adapting the approach I suggest below --in the answer to your question iii)-- in much the same way 
as is done to transform B splines into  penalized  spline. This is because the problem you try to solve can be framed as a form of (constrained!)  smoothing spline, whose objective function can then be further penalized, to force the algorithm to select the optimal knot points, in much the same way as it done in the case of P-splines. Again, see the answer to iii) below for details.

III) Is there a far more easier way to approximate any function in a
  piece-wise fashion?

A branch of statistics called non parametric boundary estimation deals with this and related problems (at least to the extend that I understand your question). The idea is to fit piece wise function to your data under the constraint that: 


*

*these piece wise function joint at pre-specified knot points,

*The fitted function is always below the real function evaluated at the knots,

*have a positive derivative (are increasing)

*(possibly) have negative second derivatives (are concave).


There are different methods depending on whether you want the piece-wise function to be linear, quadratic or cubic functions. For the last two, you can also decide whether you want to impose concavity or not (by imposing the fourth constraints as well or only the first three). Obviously, the numerical complexity of the approach depends on how many of these constraints you wish to impose but in any case it involves solving a convex optimization problem (linear program whose size depends linearly on the number of knots for the first three constraints, cone programming problem of size depending (linearly?) on the number of knot for the last constraint).  
You will find many links and R implementation at the benchmarking and npbr packages. Here is a simple example comparing DEA, FDH and LFDH (three approach in increasing order of complexity). DEA appears very close to what you are doing already, with the notable difference that it uses general trapezoid instead of box (though it can also use boxes), that in the package I link to, the code is implemented in C++ (so it will probably be a notable speed up over yours) and that the constraints I mentioned above are taken into account.
Here is a simple example to illustrate (the original CDF is in black):
library(npbr)

b0<-rexp(1000)
b1<-ecdf(b0)
x<-seq(min(b0),max(b0),l=100)
b3<-b1(x)
plot(x,b3,type="l")

z <-seq(min(x),max(x),length.out=11)
lines(z,dea_est(x,b3,z,type="dea"),lty=1,col="red")
lines(z,dea_est(x,b3,z,type="fdh"),lty=2,col="blue")
lines(z,dea_est(x,b3,z,type="lfdh"),lty=3,col="green")
legend("topleft",legend=c("dea","fdh","lfdh"),col=c("red","blue","green"),lty=1:3)


