[Robast-commits] r304 - pkg/ROptEst/inst/scripts

[noreply at r-forge.r-project.org]

Author: stamats
Date: 2009-06-13 09:37:07 +0200 (Sat, 13 Jun 2009)
New Revision: 304

Modified:
   pkg/ROptEst/inst/scripts/GammaModel.R
Log:
small shape parameters lead to functions which are hard to integrate numerically; added some code to demonstrate this ...

Modified: pkg/ROptEst/inst/scripts/GammaModel.R
===================================================================
--- pkg/ROptEst/inst/scripts/GammaModel.R	2009-06-10 13:38:10 UTC (rev 303)
+++ pkg/ROptEst/inst/scripts/GammaModel.R	2009-06-13 07:37:07 UTC (rev 304)
@@ -5,13 +5,57 @@
 options("newDevice"=TRUE)
 
 ## generates Gamma Family with 
-## scale = 2 and shape = 0.5
-G <- GammaFamily(scale = 2, shape = 0.5)
+## scale = 2 and shape = 0.1
+G <- GammaFamily(scale = 2, shape = 0.1)
 G       # show G
-plot(G) # plot of Gammad(scale = 2, shape = 0.5) and L_2 derivative
+plot(G) # plot of Gammad(scale = 2, shape = 0.1) and L_2 derivative
 distrExOptions(ErelativeTolerance = 1e-8) # increase precision for E
 checkL2deriv(G)
 
+## more precisely:
+## numerical integration gives
+E(Gammad(scale = 2, shape = 0.1), function(x) (log(x/2)-digamma(0.1))^2)
+
+## vhereas
+trigamma(0.1)
+
+## Problem is more or less caused by integration of log(x/2)^2
+## occurs for shape parameter small (<< 1)
+E(Gammad(scale = 2, shape = 0.1), function(x) log(x/2)^2)
+distrExOptions(ErelativeTolerance = 1e-10)
+E(Gammad(scale = 2, shape = 0.1), function(x) log(x/2)^2)
+distrExOptions(ErelativeTolerance = 1e-7)
+E(Gammad(scale = 2, shape = 0.1), function(x) log(x/2)^2)
+
+## result should be
+res1 <- 2*digamma(0.1)*E(Gammad(scale = 2, shape = 0.1), function(x) log(x/2))
+res2 <- digamma(0.1)^2
+trigamma(0.1) + res1 - res2
+
+fun <- function(x) log(x/2)^2*dgamma(x, scale = 2, shape = 0.1)
+integrate(fun, lower = 0, upper = Inf, rel.tol = 1e-6)
+integrate(fun, lower = 0, upper = Inf, rel.tol = 1e-3)
+integrate(fun, lower = 1e-9, upper = Inf, rel.tol = 1e-6)
+integrate(fun, lower = 1e-13, upper = Inf, rel.tol = 1e-6)
+## problem at zero!
+fun(0)
+curve(fun, from = 0, to = 1, n = 501)
+
+fun <- function(x) (log(x/2)-digamma(0.1))^2*dgamma(x, scale = 2, shape = 0.1)
+curve(fun, from = 0, to = 1, n = 501)
+
+GLIntegrate(fun, lower = 1e-9, upper = qgamma(1-1e-9, scale = 2, shape = 0.1), order = 500)
+
+
+## generates Gamma Family with 
+## scale = 2 and shape = 1
+G <- GammaFamily(scale = 2, shape = 1)
+G       # show G
+plot(G) # plot of Gammad(scale = 2, shape = 1) and L_2 derivative
+distrExOptions(ErelativeTolerance = 1e-8) # increase precision for E
+checkL2deriv(G)
+
+
 ## classical optimal IC
 IC0 <- optIC(model = G, risk = asCov())
 IC0       # show IC