#  Copyright (C) 2011-2012 Sandra Ligges
#
#  This program is free software; you can redistribute it and/or modify it under the terms
#  of the GNU General Public License as published by the Free Software Foundation;
#  either version 2 of the License, or (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
#  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
#  See the GNU General Public License for more details.
#
#  A copy of the GNU General Public License is available in the same folder as this file.


#######################################################################################################################
### Optimierungsfunktionen zur Bestimmung der Zensierungsverteilungen
#######################################################################################################################

###################################################################################################
### a) für ungleiche Zensierungsverteilungen:
###################################################################################################
### Optimierungsfunktion für die Bestimmung des Skalenparameters der Exponentialverteilung, 
### der die zu ziehenden Zensierungszeiten folgen müssen, 
### damit die Überlebenszeiten einer bestimmten Verteilung mit einer Wahrscheinlichkeit 
### von 'cens.rate' zensiert werden.
opt.cens.par.1 <- function(cens.par, cens.rate = .4, distrib = c("Exp", "Weib", "Gomp"), 
    shape = 1, scal = 1) {
    
    distrib <- match.arg(distrib)

    if((distrib == "Exp") || (distrib == "Weib")) {
        (integrate(function(x) dweibull(x, shape = shape, scale = scal^(-1/shape)) * pexp(x, cens.par), 0, Inf)$value
         - cens.rate)^2
    }
    else {
        (integrate(function(x) dgompertz(x, alpha = shape, theta = scal) * pexp(x, cens.par), 0, 1000)$value
        - cens.rate)^2
    }
}



###################################################################################################
### b) für gleiche Zensierungsverteilungen:
###################################################################################################
### Optimierungsfunktion für die Bestimmung des Skalenparameters der Exponentialverteilung, 
### der die zu ziehenden Zensierungszeiten folgen müssen, 
### damit der gruppenübergreifende Zensierungsanteil unter den Überlebenszeiten zweier Gruppen 
### mit unterschiedlicher Verteilung bei 'cens.rate' liegt.
opt.cens.par.2 <- function(cens.par, cens.rate = .4, distrib = c("Exp", "Weib", "Gomp"), 
    shape = c(1, 1), scal = c(1, 1), w1 = 0.5, w2 = 0.5) {
    
    distrib <- match.arg(distrib)
    shape1 <- shape[1]; shape2 <- shape[2]
    scal1 <- scal[1]; scal2 <- scal[2]

    if((distrib == "Exp") || (distrib == "Weib")) {
        (w1 * integrate(function(x) dweibull(x, shape = shape1, scale = scal1^(-1/shape1)) * 
            pexp(x, cens.par), 0, Inf)$value
        + w2 * integrate(function(x) dweibull(x, shape = shape2, scale = scal2^(-1/shape2)) * 
            pexp(x, cens.par), 0, Inf)$value
        - cens.rate)^2
    }
    else {
        (w1 * integrate(function(x) dgompertz(x, alpha = shape1, theta = scal1) * pexp(x, cens.par), 0, 500)$value
        + w2 * integrate(function(x) dgompertz(x, alpha = shape2, theta = scal2) * pexp(x, cens.par), 0, 500)$value
        - cens.rate)^2
    }
}