关于矩阵求逆奇异值的问题

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

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

立即领取

感谢您参与论坛问题回答

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

+2 论坛币

1. #NO.5
   
2. data <- read.table(file.choose(),header=TRUE)
   
3. names(data)
   
4. x <- data$probability
   
5. # a) Estimate by Moments
   
6. x_mean <- mean(x) ;x_var <- var(x);x_mean;x_var
   
7. alpha0 <- (x_mean*x_var)/(x_mean^2*(1-x_mean)-x_var);alpha0
   
8. beta0 <- alpha0*(1-x_mean)/x_mean;beta0
   
9. 
   
10. # b) Estimate by Likelihood function
    
11. #定义对数似然函数
    
12. n <- nrow(as.matrix(x));n
    
13. alpha_old <- alpha0; beta_old <- beta0
    
14. LL_func <- function(x, alpha=alpha_old, beta=bata_old){
    
15.   LL <- -n*log(beta(alpha_old,beta_old))+(alpha_old-1)*sum(x)+(beta_old-1)*sum(1-x)
    
16. }
    
17. #用Newton-Raphson数值解的方法，求最大似然函数的解。
    
18. tol0 <- 1e-3;kmax <- 500; k = 1;
    
19. tol <-1
    
20. while(tol > tol0 | k <= kmax){
    
21.   delta_alpha <- 0.1*alpha_old;delta_beta <- 0.1*beta_old
    
22.   par_old <- matrix(alpha_old,beta_old,ncol=1)
    
23.   f1 <- LL_func(x,alpha_old+delta_alpha,beta_old+delta_beta)
    
24.   f2 <- LL_func(x,alpha_old+delta_alpha,beta_old-delta_beta)
    
25.   f3 <- LL_func(x,alpha_old-delta_alpha,beta_old+delta_beta)
    
26.   f4 <- LL_func(x,alpha_old-delta_alpha,beta_old-delta_beta)
    
27.   f5 <- LL_func(x,alpha_old,beta_old-delta_beta)
    
28.   f6 <- LL_func(x,alpha_old+delta_alpha,beta_old)
    
29.   #定义梯度阵和Hasen矩阵
    
30.   G <- matrix((f1+f2-f3-f4)/(4*delta_alpha),(f1+f3-f2-f4)/(4*delta_beta),ncol = 1)
    
31.   H <- matrix(c((f1+f3-2*f5)/(2*(delta_alpha^2)),
    
32.               (f1+f4-f2-f3)/(4*delta_alpha*delta_beta),
    
33.               (f1+f4-f2-f3)/(4*delta_alpha*delta_beta),
    
34.               (f1+f3-2*f6)/(2*(delta_beta^2))),nrow=2,ncol=2)
    
35.   par_new <-par_old - solve(H) %*% G
    
36.   tol <- t(par_old-par_new) %*% (par_old-par_new)
    
37.   par_old <- par_new
    
38.   alpha_old <- par_old[1,1];beta_old <- par_old[2,1]
    
39.   delta_alpha <- 0.1*alpha_old;delta_beta <- 0.1*beta_old
    
40.   k=k+1
    
41.   print (par_new)
    
42. }

复制代码

influenza.txt (12.16 KB)

2013-11-7 16:52:11 上传

数据集

坐等神人。

扫码加我 拉你入群

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

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

分享0 收藏0 回帖

关键词：奇异值 Data

Error in solve.default(H) :
  Lapack routine dgesv: system is exactly singular: U[1,1] = 0
> H
     [,1] [,2]
[1,]    0    0
[2,]    0    0
我觉得可能是H值得问题，H不能够全是0
然后，我有看下
> f1
[1] 423.304
> f2
[1] 423.304
> f3
[1] 423.304
> f4
[1] 423.304
> f5
[1] 423.304
> f6
[1] 423.304

这几个值都相等，符合你的要求吗，如果不符合，可能是你自定义的LL_func函数有问题