在R中解方程:

huangapple go评论114阅读模式
英文:

Solving an equation in r

问题

假设我有以下方程:

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

我有所有变量的值,除了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的值?

我尝试使用unirootoptimise函数,但仍然没有得到答案。也许我使用了错误的函数?

请尝试使用uniroot函数来解决这个方程以获得sigma的值。uniroot函数是用于单变量非线性方程求解的强大工具。以下是如何使用它的示例:

# 定义函数
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

# 使用uniroot求解
result <- uniroot(sigmagm, interval = c(0, 1)) # 调整区间以适应sigma的可能值范围

# 结果
sigma_value <- result$root
cat("Sigma的值为:", sigma_value)

这将尝试找到满足sigmagm函数的根,即sigma的值。记得根据问题的背景调整区间,以包含可能的sigma值范围。

英文:

Suppose I have the following equation:

sigmagm &lt;- 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?

答案1

得分: 3

为了获得反函数,我们首先需要获得 sigmagm 函数的定义域。可以得知 sigmagm(0,exp(-7.575783407)] 内定义,该定义域之外的任何值都会产生 NA。同时注意到在区间 (0, exp(-31.7064)) 内该函数有一个根,而在大于 exp(-31.7064) 的区间内该函数有两个根。

在进行反函数时,我们要注意到计算机在计算非常小的数字时存在限制。因此,该函数的取值范围为 (18.30267, 162),综合考虑所有这些因素,反函数可以计算如下:

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

注意,对于任何小于最小值的 x 值,我们会返回最小值,而对于任何 x>162,我们会将 x 设为 162

因此:

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))
  
}
英文:

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 &lt;- function(x){
  upper_limit &lt;- exp(-7.575783407)
  fn_upper &lt;-  sigmagm(upper_limit)
  div &lt;- sigmagm(1.698581080212673947126e-14)
  fn &lt;- function(i, x){
    if(i[1]&gt;upper_limit || i[1]&lt;=0) NA
    else abs(sigmagm(i[1]) - x)
  }
  suppressWarnings(c(optim(2e-30, fn, x = x)[[1]],
    if(x &lt;= 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&gt;162 we set x<-162.

Thus:

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

sigmagm &lt;- 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))
  
}

huangapple
  • 本文由 发表于 2023年7月7日 03:10:47
  • 转载请务必保留本文链接:https://go.coder-hub.com/76631893.html
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定