## distribucio'n del beneficio
dis <- function (C = 6,             # encargos
                 Dmin = 4,          # demanda mi'nima
                 N = 4,             # demanda = Dmin +
                 p = 3/4,           # + Binomial (N; p)
                 cos = 4,           # costo unitario
                 pvp = 5)           # precio de venta
{                                   
    D <- Dmin : (Dmin+N)            # valores posibles de la demanda
    P <- dbinom(0:N, N, p)          # probabilidades de la demanda
    V <- pmin(D, C)                 # cantidades vendidas
    B <- pvp * V - cos * C          # beneficio
    P <- by(P, B, sum)              # probabilidades del beneficio
    B <- as.numeric (names (P))     # modalidades del beneficio
    list(val=B, prob=as.numeric(P)) # distribucio'n: valores y probabilidades
}

mediaB <- function(...)                    # beneficio esperado
{
    distro <- dis(...)                     # distribucio'n de probabilidades
    weighted.mean(distro$val, distro$prob) # media ponderada
}

varB <- function(...)                      # varianza del beneficio
{
    distro <- dis(...)
    weighted.mean(distro$val^2, distro$prob) - mediaB(...)^2
}

mejor <- function (Dmin=4, N=4, ...)       # nu'mero o'ptimo de encargos
{
    encargables <- Dmin : (Dmin+N)         # valores razonables 
    beneficios  <- sapply(encargables,     # beneficios esperados
                          function(C) mediaB(C, Dmin, N, ...))
    names(beneficios) <- encargables       # para identificar el C o'ptimo
    beneficios [which.max (beneficios)]    # o'ptimo y su C asociado
}

plot(xx <- 4:24, sapply(xx, function(C) mediaB(N=20,C=C)))
mejor(N=20)