CODEWITHSUNDEEP
requationequation-solving
Farrelyto Theodorus
Farrelyto Theodorus

Reputation: 21

Solving an equation in r

Suppose I have the following equation:

sigmagm <- function(sigma){
  
(-sigma*t) - (alpha_h*t) - ((beta1_dot_h/beta2_h) * (exp(beta2_h*(x+t)) - exp(beta2_h*x))) - ((beta1_dot_h/beta2_h) * exp(beta2_h*x)*(1 - (1+gama) * exp(beta2_h*t))) + (alpha_dd + (1+gama) * alpha_h) * t + ((beta1_dot_dd/beta2_dd) * exp(beta2_dd * (x+t))) - (gama * (beta1_dot_h/beta2_h) * exp(beta2_h * (x+t)) * (1-(beta2_h * (t/2)))) - ((beta1_dot_dd/beta2_dd) * exp(beta2_dd * (x+(t/2))) * (1 - beta2_dd*(t/2))) - log(sigma) - log(exp((gama*alpha_h) + alpha_dd - sigma + (gama*beta1_dot_h*exp(beta2_h * (x+(t/2)))) + (beta1_dot_dd * exp(beta2_dd * (x+(t/2)))*t)) - 1) - log((gama*alpha_h) + alpha_dd - sigma + (gama * beta1_dot_h * exp(beta2_h * (x+(t/2)))) + (beta1_dot_dd*exp(beta2_dd * (x+(t/2))))) - log((1 - prev)/prev)
  
  return(sigmagm)
}

and I have the value of all variables except for the sigma

alpha_h = -0.001188
beta1_h = -5.85382
beta2_h = 0.43633
alpha_dd = -0.002006
beta1_dd = -6.24563
beta2_dd = 0.01495
beta1_dot_dd = exp(beta1_dd)
beta1_dot_h = exp(beta1_h)
prev = 0.057115867
t = 5
x = 15
gama = 0

How can I solve the equation to obtain the value of sigma?

I've been trying using uniroot and optimise function but I still don't get the answer. Maybe I use the wrong function?

Upvotes: 0

Views: 61

Answers (1)

Onyambu
Onyambu

Reputation: 79328

To obtain the inverse function, we first need to obtain the domain of the sigmagm function. it follows that the sigmagm is defined within: (0,exp(-7.575783407)] Anything ouside of this domain produces NA. Notice also that between (0, exp(-31.7064)) the function has one root while anything above exp(-31.7064) the function has two roots. When doing the inversion, we note that the computer is limited with the computation of very small numbers. The range thus for the function is (18.30267, 162) Taking all these into account, the inverse can be computed as:

Inv_Sigmagm <- function(x){
  upper_limit <- exp(-7.575783407)
  fn_upper <-  sigmagm(upper_limit)
  div <- sigmagm(1.698581080212673947126e-14)
  fn <- function(i, x){
    if(i[1]>upper_limit || i[1]<=0) NA
    else abs(sigmagm(i[1]) - x)
  }
  suppressWarnings(c(optim(2e-30, fn, x = x)[[1]],
    if(x <= div) optim(upper_limit, fn, x = x)[[1]]))
  
}

Notice that for any x value below tha minimum, we will return the minimum while for any x>162 we set x<-162.

Thus:

Inv_Sigmagm(50)
[1] 2.146653e-18
sigmagm(Inv_Sigmagm(50))
[1] 50

sigmagm <- function(sigma){
  alpha_h = -0.001188
  beta1_h = -5.85382
  beta2_h = 0.43633
  alpha_dd = -0.002006
  beta1_dd = -6.24563
  beta2_dd = 0.01495
  beta1_dot_dd = exp(beta1_dd)
  beta1_dot_h = exp(beta1_h)
  prev = 0.057115867
  t = 5
  x = 15
  gama = 0
  
  ((-sigma*t) - (alpha_h*t) - 
    ((beta1_dot_h/beta2_h) * (exp(beta2_h*(x+t)) - exp(beta2_h*x))) - 
    ((beta1_dot_h/beta2_h) * exp(beta2_h*x)*(1 - (1+gama) * exp(beta2_h*t))) + 
    (alpha_dd + (1+gama) * alpha_h) * t + 
    ((beta1_dot_dd/beta2_dd) * exp(beta2_dd * (x+t))) - 
    (gama * (beta1_dot_h/beta2_h) * exp(beta2_h * (x+t)) * (1-(beta2_h * (t/2)))) - 
    ((beta1_dot_dd/beta2_dd) * exp(beta2_dd * (x+(t/2))) * (1 - beta2_dd*(t/2))) - 
    log(sigma) - 
    log(exp((gama*alpha_h) + alpha_dd - sigma +
              (gama*beta1_dot_h*exp(beta2_h * (x+(t/2)))) + 
              (beta1_dot_dd * exp(beta2_dd * (x+(t/2)))*t)) - 1) -
    log((gama*alpha_h) + alpha_dd - sigma + 
          (gama * beta1_dot_h * exp(beta2_h * (x+(t/2)))) +
          (beta1_dot_dd*exp(beta2_dd * (x+(t/2))))) - 
    log((1 - prev)/prev))
  
}

Upvotes: 3

Related Questions