R/calcul.bi.R

Defines functions calcul.bi

#' @import ggplot2
#' @import rstan
#' @import methods
#' @import stats
#' @import graphics
#' @import grDevices
#' @import stats4
calcul.bi <-
function(wm, meanbeta, a = NULL, model, tau){
    xistar <- calcul.xi(wm, meanbeta, a, model)
    if(model == "logistic"){
        f <- function(b, xj, xjbis){1/(1+exp(-a-exp(b)*xj)) + 1/(1+exp(-a-exp(b)*xjbis)) - 2*tau}
        lesb <- c()
        for(i in 1:(length(xistar)-1)){
            b <- uniroot(f, c(-1000000, 100000), xj = xistar[i], xjbis = xistar[i + 1] )$root
            lesb <- c(lesb, b)
        }
    }else{
        if(model == "power_log"){
            f <- function(b, xj, xjbis){xj^exp(b) + xjbis^exp(b) - 2*tau}
            lesb <- c()
            for(i in 1:(length(xistar) - 1)){
                b <- uniroot(f, c(-1000, 1000), xj = xistar[i], xjbis = xistar[i + 1] )$root
                lesb <- c(lesb, b)
            }
            bl <- log(log(tau + 0.05)/log(wm[1]))
            bu <- log(log(tau - 0.05)/log(wm[length(wm)]))
            lesb <- c(bl, lesb, bu)
        }else{
            if(model == "power"){
                f <- function(b, xj, xjbis){xj^b + xjbis^b - 2*tau}
                lesb <- c()
                for(i in 1:(length(xistar) - 1)){
                    b <- uniroot(f, c(-1000, 1000), xj = xistar[i], xjbis = xistar[i + 1] )$root
                    lesb <- c(lesb, b)
                }
                bl <- log(tau + 0.05)/log(wm[1])
                bu <- log(tau - 0.05)/log(wm[length(wm)])
                lesb <- c(bl, lesb, bu)
            }else{stop("Error: Please insert a valid model")}
        }
    }
    return(lesb)
}

Try the dfped package in your browser

Any scripts or data that you put into this service are public.

dfped documentation built on May 2, 2019, 8:36 a.m.