怎么用R求解一个方程组。。。

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

这个是方程组，求解的是均值和标准差
我的程序是这样的
x<-rnorm(100,0,1)
n<-length(x)
funs<-function(x){
f<-c(sum(x-x[1])/x[2]^2,(-n/2*x[2]^2)+sum((x-x[1])^2)/2*x[2]^4)
J<-matrix(c(-n/x[2]^2,-2*sum(x-x[1])/x[2]^3,-sum(x-x[1])/x[2]^4,(n/x[2]^3)-2*sum((x-x[1])^2)/x[2]^5),nrow=2,byrow=T)
list(f=f,J=J)
}
Newtons<-function(fun,x,ep=1e-5,it_max=2000){
index<-0;k<-1
while(k<=it_max){
x1<-x;obj<-fun(x);x<-x-solve(obj$J,obj$f);
norm<-sqrt((x-x1)%*%(x-x1))
if(norm<ep){
index<-1;break
}
k<-k+1
}
obj<-fun(x);
list(root=f,it=k,index=index,FunVal=obj$f)
}
Newtons(funs,c(1,1))


但是每次都出现这个结果Error in solve.default(obj$J, obj$f) :
  system is computationally singular: reciprocal condition number = 6.3691e-34
求大神解答是哪里出问题了？

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏2 回帖

关键词：方程组 标准差 用R求解方程组

我把代码写的比较原始，改动的地方有两个：
1.在函数funs里把x[1]和x[2]改成了z[1],z[2]，同时将sigma平方做为未知数
2.在函数newtons里把norm改动了一下

1. > rm(list=ls())
   
2. > ls()
   
3. character(0)
   
4. > set.seed(91)
   
5. > x<-rnorm(100,0,1)
   
6. > n<-length(x)
   
7. > funs<-function(z){
   
8. +     J <- matrix(1:4,nrow=2)
   
9. +     f <- matrix(1:2,nrow=2)
   
10. +     f[1,1]<-sum(x-z[1])/z[2]
    
11. +     f[2,1] <- -n/2*z[2]+sum((x-z[1])^2)/(2*z[2]^2)
    
12. +
    
13. +     J[1,1] <- -n/z[2]
    
14. +     J[2,1] <- -sum(x-z[1])/z[2]^2
    
15. +     J[1,2] <- -sum(x-z[1])/z[2]^2
    
16. +     J[2,2] <- n/(2*z[2]^2)-sum((x-z[1])^2)/z[2]^3
    
17. +     list(f=f,J=J)
    
18. + }
    
19. > Newtons<-function(fun,z,ep=0.00001,it_max=2000){
    
20. +     index<-0;k<-1
    
21. +     while(k<=it_max){
    
22. +         z1<-z;obj<-fun(z);z<-z-solve(obj$J,obj$f)
    
23. +         norm<-(abs(z1[1]-z[1])<ep)&(abs(z1[2]-z[2])<ep)
    
24. +         #stopifnot(any(is.na(norm)))
    
25. +         #print(obj)
    
26. +         #cat("root",z)
    
27. +         if(norm){
    
28. +             index<-1
    
29. +             stop
    
30. +         }
    
31. +         k<-k+1
    
32. +     }
    
33. +     #obj<-fun(z);
    
34. +     list(root=z,it=k,index=index,FunVal=obj$f)
    
35. + }
    
36. > Newtons(funs,c(1,1))
    
37. $root
    
38.             [,1]
    
39. [1,] -0.08441826
    
40. [2,]  1.17773854
    
41. 
    
42. $it
    
43. [1] 2001
    
44. 
    
45. $index
    
46. [1] 0
    
47. 
    
48. $FunVal
    
49.               [,1]
    
50. [1,] -2.691703e-15
    
51. [2,]  5.394353e+01
    
52. 
    
53. > mean(x);var(x)
    
54. [1] -0.08441826
    
55. [1] 1.348797

复制代码

均值的未知数能得到真实值，但是方差未知数差距很大
迭代一百万次后的结果

1. > Newtons(funs,c(1,1),it_max=1000000)
   
2. $root
   
3.             [,1]
   
4. [1,] -0.08441826
   
5. [2,]  1.17773854
   
6. 
   
7. $it
   
8. [1] 1000001
   
9. 
   
10. $index
    
11. [1] 0
    
12. 
    
13. $FunVal
    
14.               [,1]
    
15. [1,] -2.691703e-15
    
16. [2,]  5.394353e+01

复制代码

他们是第一个方程的解，但不是第二个方程的解，估计是设置方程的时候敲错了符号或者其他什么原因

找到错误了，是公式书写错误，改正后的结果

1. > rm(list=ls())
   
2. > ls()
   
3. character(0)
   
4. > set.seed(91)
   
5. > x<-rnorm(100,0,1)
   
6. > n<-length(x)
   
7. > funs<-function(z){
   
8. +     J <- matrix(1:4,nrow=2)
   
9. +     f <- matrix(1:2,nrow=2)
   
10. +     f[1,1]<-sum(x-z[1])/z[2]
    
11. +     f[2,1] <- -n/(2*z[2])+sum((x-z[1])^2)/(2*z[2]^2)
    
12. +
    
13. +     J[1,1] <- -n/z[2]
    
14. +     J[2,1] <- -sum(x-z[1])/z[2]^2
    
15. +     J[1,2] <- -sum(x-z[1])/z[2]^2
    
16. +     J[2,2] <- n/(2*z[2]^2)-sum((x-z[1])^2)/z[2]^3
    
17. +     list(f=f,J=J)
    
18. + }
    
19. > Newtons<-function(fun,z,ep=0.00001,it_max=2000){
    
20. +     index<-0;k<-1
    
21. +     while(k<=it_max){
    
22. +         z1<-z
    
23. +         obj<-fun(z)
    
24. +         z<-z-solve(obj$J,obj$f)
    
25. +         norm<-(abs(z1[1]-z[1])<ep)&(abs(z1[2]-z[2])<ep)
    
26. +         #norm <- sum(z1[2]-z[1])<ep
    
27. +         #stopifnot(any(is.na(norm)))
    
28. +         #print(obj)
    
29. +         #cat("root",z)
    
30. +         if(norm){
    
31. +             index<-1
    
32. +             break
    
33. +         }
    
34. +         k<-k+1
    
35. +     }
    
36. +     #obj<-fun(z);
    
37. +     list(root=z,it=k,index=index,FunVal=obj$f)
    
38. + }
    
39. > Newtons(funs,c(0.9,0.96),it_max=10000)
    
40. $root
    
41.             [,1]
    
42. [1,] -0.08441826
    
43. [2,]  1.33530886
    
44. 
    
45. $it
    
46. [1] 8
    
47. 
    
48. $index
    
49. [1] 1
    
50. 
    
51. $FunVal
    
52.             [,1]
    
53. [1,] 1.27524e-07
    
54. [2,] 8.48539e-06
    
55. 
    
56. > mean(x);sd(x)
    
57. [1] -0.08441826
    
58. [1] 1.161377
    
59. > var(x)
    
60. [1] 1.348797
    
61. >

复制代码

这个可以用了，嘿嘿

nuomin 发表于 2014-3-6 21:58
我把代码写的比较原始，改动的地方有两个：
1.在函数funs里把x[1]和x[2]改成了z[1],z[2]，同时将sigma平方 ...

太感谢您了。。。

谢谢分享