- #NO.5
- data <- read.table(file.choose(),header=TRUE)
- names(data)
- x <- data$probability
- # a) Estimate by Moments
- x_mean <- mean(x) ;x_var <- var(x);x_mean;x_var
- alpha0 <- (x_mean*x_var)/(x_mean^2*(1-x_mean)-x_var);alpha0
- beta0 <- alpha0*(1-x_mean)/x_mean;beta0
- # b) Estimate by Likelihood function
- #定义对数似然函数
- n <- nrow(as.matrix(x));n
- alpha_old <- alpha0; beta_old <- beta0
- LL_func <- function(x, alpha=alpha_old, beta=bata_old){
- LL <- -n*log(beta(alpha_old,beta_old))+(alpha_old-1)*sum(x)+(beta_old-1)*sum(1-x)
- }
- #用Newton-Raphson数值解的方法，求最大似然函数的解。
- tol0 <- 1e-3;kmax <- 500; k = 1;
- tol <-1
- while(tol > tol0 | k <= kmax){
- delta_alpha <- 0.1*alpha_old;delta_beta <- 0.1*beta_old
- par_old <- matrix(alpha_old,beta_old,ncol=1)
- f1 <- LL_func(x,alpha_old+delta_alpha,beta_old+delta_beta)
- f2 <- LL_func(x,alpha_old+delta_alpha,beta_old-delta_beta)
- f3 <- LL_func(x,alpha_old-delta_alpha,beta_old+delta_beta)
- f4 <- LL_func(x,alpha_old-delta_alpha,beta_old-delta_beta)
- f5 <- LL_func(x,alpha_old,beta_old-delta_beta)
- f6 <- LL_func(x,alpha_old+delta_alpha,beta_old)
- #定义梯度阵和Hasen矩阵
- G <- matrix((f1+f2-f3-f4)/(4*delta_alpha),(f1+f3-f2-f4)/(4*delta_beta),ncol = 1)
- H <- matrix(c((f1+f3-2*f5)/(2*(delta_alpha^2)),
- (f1+f4-f2-f3)/(4*delta_alpha*delta_beta),
- (f1+f4-f2-f3)/(4*delta_alpha*delta_beta),
- (f1+f3-2*f6)/(2*(delta_beta^2))),nrow=2,ncol=2)
- par_new <-par_old - solve(H) %*% G
- tol <- t(par_old-par_new) %*% (par_old-par_new)
- par_old <- par_new
- alpha_old <- par_old[1,1];beta_old <- par_old[2,1]
- delta_alpha <- 0.1*alpha_old;delta_beta <- 0.1*beta_old
- k=k+1
- print (par_new)
- }
坐等神人。 |
熟悉论坛请点击新手指南
京公网安备 11010802022788号
论坛法律顾问：王进律师
知识产权保护声明
免责及隐私声明
