Solving an equation in r

Question

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?

If the equation above calculates ```sigmagm``` from these variables and ```sigma```, then is ```sigmagm``` given in this case? — Archeologist, Jul 06 '23 at 19:24
`sigmagm` you provided is a function that takes a single parameter `sigma`. There is no equation to be solved here! Other variables used in the expression (not equation) are not defined within the function itself but are assumed to be defined elsewhere in the code or in the global environment. — Byte Ninja, Jul 06 '23 at 19:27
Are you sure that function has roots? Also, you should remove the `return` instruction, all of it, for the function to make sense, like this you are returning its argument. — Rui Barradas, Jul 06 '23 at 19:36
first you only solve an equation: eg `3x+1=4` solve for `x`. you note that you have the function value `4` and one missing parameter. in youe case, you have not provided a function value to be equated unto. how do you want to proceed? — Onyambu, Jul 06 '23 at 20:44
The function appears not to have any roots. It does have a minimum value of ~18.30628 at `sigma = 0.0002687089`. Try `curve(sigmagm, 1e-6, 5e-4)` and `optimize(sigmagm, c(1e-6, 5e-4))`. — jblood94, Jul 06 '23 at 20:57

score 3 · Answer 1 · answered Jul 06 '23 at 22:04

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))
  
}

Solving an equation in r

1 Answers1