提交: d98741928ecb2febd67065825ddfa9d0c8d43ca1 - myprojects-hg-reborn (hg)

Commit MetaInfo

• 修订版: d98741928ecb2febd67065825ddfa9d0c8d43ca1 (tree)
• 时间: 2017-12-21 19:57:18
• 作者: Lorenzo Isella <lorenzo.isella@gmai...>
• Commiter: Lorenzo Isella

Log Message

A code showing how to fit a beta distributions using a couple of methods and also explicitly writing the log-likelyhood function.

更改概述

• add: R-codes/beta_fitting_test.R

差异

diff -r e0b6162efbb0 -r d98741928ecb R-codes/beta_fitting_test.R
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/R-codes/beta_fitting_test.R	Thu Dec 21 11:57:18 2017 +0100
@@ -0,0 +1,87 @@
+rm(list=ls())
+
+library(MASS)
+library(stats4)
+
+
+dbeta1 <- function(x, shape, ...)
+       dbeta(x, shape, shape, ...)
+
+
+dbeta2 <- function(x, shape){
+
+
+    res <- x^(shape-1)*(1-x)^(shape-1)/beta(shape, shape)
+
+    return(res)
+
+}
+
+
+
+LL1 <- function(shape){
+
+
+    R <- dbeta1(x, shape)
+
+
+    res <- -sum(log(R))
+
+    return(res)
+
+
+}
+
+
+
+
+LL2 <- function(shape){
+
+
+    R <- dbeta2(x, shape)
+
+
+    res <- -sum(log(R))
+
+    return(res)
+
+
+}
+
+
+
+
+set.seed(124)
+
+x <-rbeta(1000, 0.2, 0.2)
+
+    
+fit_dbeta1 <- fitdistr( x , dbeta1, start=list(shape=0.5) , method="Brent", lower=c(0), upper=c(1))
+
+print("estimate of shape from fit_dbeta1 is")
+print(fit_dbeta1$estimate)
+
+
+
+
+fit_dbeta2 <- fitdistr( x , dbeta2, start=list(shape=0.5) , method="Brent", lower=c(0), upper=c(1))
+
+print("estimate of shape from fit_dbeta2 is")
+print(fit_dbeta2$estimate)
+
+
+
+fit_LL1 <- mle(LL1, start=list(shape=0.5))
+
+print("estimate of from fit_LL1")
+print(coef(fit_LL1))
+
+
+## this does not work
+fit_LL2 <- mle(LL2, start=list(shape=0.5))
+print(coef(fit_LL2))
+
+
+
+
+print("So far so good")