英文:
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的值?
我尝试使用uniroot
和optimise
函数,但仍然没有得到答案。也许我使用了错误的函数?
请尝试使用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 <- 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 <- 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))
}
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库,让每个人都能够通过互相帮助和分享经验来进步。
评论