求解一个简单的非线性系统
Solving a simply non linear system
例如我有这个数据:
x-c(73,6,77,81,91,120,150,61,65,68,18,20,23,12,14,18,23,26
+26,27,2,3,3,40,41,41,6,10,11,12,37,38,38,6,73,6,51)
我想计算伽玛分布的 a=shape 和 b=scale 参数。我想解决这个非线性系统
a*b=m1
a*b^2+(a^2)*(b^2)=m2
m1 和 m2 是这些:
m1<-sum(x)/length(x)
m2<-sum((x)^2)/length(x)
我可以用手和计算器解决它,但我想知道如何用 R 立即解决这个问题
R 可以轻松做到这一点
library(nleqslv)
f <- function(x) {
a<-x[1]
b<-x[2]
c(a*b-m1,a*b^2+(a^2)*(b^2)-m2)
}
nleqslv(c(1,30), f)
输出应如下所示:
$`x`
[1] 1.286486 30.595840
$fvec
[1] -9.663381e-13 -1.396074e-10
$termcd
[1] 1
$message
[1] "Function criterion near zero"
$scalex
[1] 1 1
$nfcnt
[1] 11
$njcnt
[1] 2
$iter
[1] 10
您可以通过提供渐变使事情变得更稳健。当然,R 也可以直接估计伽马分布的参数(例如来自 MASS 包的 fitdistr
)。
例如我有这个数据:
x-c(73,6,77,81,91,120,150,61,65,68,18,20,23,12,14,18,23,26
+26,27,2,3,3,40,41,41,6,10,11,12,37,38,38,6,73,6,51)
我想计算伽玛分布的 a=shape 和 b=scale 参数。我想解决这个非线性系统
a*b=m1
a*b^2+(a^2)*(b^2)=m2
m1 和 m2 是这些:
m1<-sum(x)/length(x)
m2<-sum((x)^2)/length(x)
我可以用手和计算器解决它,但我想知道如何用 R 立即解决这个问题
R 可以轻松做到这一点
library(nleqslv)
f <- function(x) {
a<-x[1]
b<-x[2]
c(a*b-m1,a*b^2+(a^2)*(b^2)-m2)
}
nleqslv(c(1,30), f)
输出应如下所示:
$`x`
[1] 1.286486 30.595840
$fvec
[1] -9.663381e-13 -1.396074e-10
$termcd
[1] 1
$message
[1] "Function criterion near zero"
$scalex
[1] 1 1
$nfcnt
[1] 11
$njcnt
[1] 2
$iter
[1] 10
您可以通过提供渐变使事情变得更稳健。当然,R 也可以直接估计伽马分布的参数(例如来自 MASS 包的 fitdistr
)。