[Distr-commits] r1243 - in branches/distr-2.8/pkg/distrMod: . R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Aug 5 17:55:36 CEST 2018
Author: ruckdeschel
Date: 2018-08-05 17:55:36 +0200 (Sun, 05 Aug 2018)
New Revision: 1243
Added:
branches/distr-2.8/pkg/distrMod/R/asCvMVarianceQtl.R
Modified:
branches/distr-2.8/pkg/distrMod/NAMESPACE
branches/distr-2.8/pkg/distrMod/R/0distrModUtils.R
branches/distr-2.8/pkg/distrMod/R/MCEstimator.R
branches/distr-2.8/pkg/distrMod/R/MDEstimator.R
branches/distr-2.8/pkg/distrMod/R/MLEstimator.R
branches/distr-2.8/pkg/distrMod/R/SimpleL2ParamFamilies.R
branches/distr-2.8/pkg/distrMod/R/internalMleCalc.R
branches/distr-2.8/pkg/distrMod/inst/NEWS
branches/distr-2.8/pkg/distrMod/man/MDEstimator.Rd
branches/distr-2.8/pkg/distrMod/man/internalmleHelpers.Rd
branches/distr-2.8/pkg/distrMod/man/internals.Rd
Log:
[distrMod] branch 2.8:
+ distinguish two cases for CvMMDEstimator: mu = emp. cdf (default) and mu = current best fit model distribution
(controlled by argument muDatOrMod = c("Dat","Mod")) => consistency between estimate and asyCov
+ added some theory/references to help file to MD estimators
bug fixes
+ discovered some issues with local variables in L2Families (global values were used instead...)
+ default for asCov of CvMMDEstimator was inconsistent: the estimator was using emp.cdf, but
the asyCov was using mu the current best fit model distribution
+ in the wrappers to MDEstimator: CvMMDEstimator, KolmogorovMDEstimator, TotalVarMDEstimator,
HellingerMDEstimator, we had the "wrong" call in slot estimate.call
+ in the last commit, forgot to change param to param.0 when calling L2ParamFamily in
SimpleL2ParamFamilies.R (when res <- L2ParamFamily(....))
+ in the last commit, mixed prob.0 and prob.1 in line 33
under the hood
+ replaced integration for AbscontDistribution(s) in .CvMMDCovariance by integration on quantile scale
=> CvMMDEstimator now works with variances even for Gamma distributions for shape < 1 ...
+ .process.meCalcRes gains arg "x" to be able to pass on emp.CDF for mu in CvMMDEstimator
if arg asvar.fct of MCEstimator has "x" in formals the observations x are passed on to asvar.fct,
otherwise they are not; correspondingly "x" is passed on to .process.meCalcRes in
MCEstimator(), MDEstimator(), MLEstimator().
+ old .CvMMDCovariance() becomes .oldCvMMDCovariance
+ new wrapper .CvMMDCovarianceWithMux which uses emp cdf as mu
+ new wrappe CvMDist2 which by default uses model distribution as mu
+ CvMMDEstimator gains argument muDatOrMod = c("Dat","Mod") to distinguish two cases
+ moved code to .[old]CvMMDCovariance from 0distrModUtils.R to new file asCvMVarianceQtl.R
Modified: branches/distr-2.8/pkg/distrMod/NAMESPACE
===================================================================
--- branches/distr-2.8/pkg/distrMod/NAMESPACE 2018-08-05 15:52:18 UTC (rev 1242)
+++ branches/distr-2.8/pkg/distrMod/NAMESPACE 2018-08-05 15:55:36 UTC (rev 1243)
@@ -94,4 +94,5 @@
export("L2LocationUnknownScaleFamily", "L2ScaleUnknownLocationFamily")
export("meRes", "get.criterion.fct")
export("addAlphTrsp2col")
-export(".deleteDim",".isUnitMatrix",".CvMMDCovariance")
+export(".deleteDim",".isUnitMatrix",
+ ".CvMMDCovariance", ".oldCvMMDCovariance", ".CvMMDCovarianceWithMux")
Modified: branches/distr-2.8/pkg/distrMod/R/0distrModUtils.R
===================================================================
--- branches/distr-2.8/pkg/distrMod/R/0distrModUtils.R 2018-08-05 15:52:18 UTC (rev 1242)
+++ branches/distr-2.8/pkg/distrMod/R/0distrModUtils.R 2018-08-05 15:55:36 UTC (rev 1243)
@@ -90,318 +90,6 @@
}
-.CvMMDCovariance<- function(L2Fam, param, mu = distribution(L2Fam),
- withplot = FALSE, withpreIC = FALSE,
- N = getdistrOption("DefaultNrGridPoints")+1,
- rel.tol=.Machine$double.eps^0.3,
- TruncQuantile = getdistrOption("TruncQuantile"),
- IQR.fac = 15,
- ...){
-
- # preparations:
-
- N1 <- 2*N+1
- odd <- (1:N1)%%2==1
-
- param0 <- L2Fam at param
- dim0 <- dimension(param0)
-# print(param0)
- paramP <- param0
- paramP at main <- main(param)
- paramP at trafo <- diag(dim0)
-# print(paramP)
- L2Fam <- modifyModel(L2Fam, paramP)
-
-# print(L2deriv(L2Fam)[[1]]@Map)
- distr <- L2Fam at distribution
-
- ### get a sensible integration range:
- low0 <- q.l(distr)(TruncQuantile)
- up0 <- q.l(distr)(TruncQuantile, lower.tail = FALSE)
- m0 <- median(distr); s0 <- IQR(distr)
- low1 <- m0 - IQR.fac * s0
- up1 <- m0 + IQR.fac * s0
- low <- max(low0,low1); up <- min(up0,up1)
-
- ### get a sensible integration range:
- if(missing(mu)) mu <- distr
- low0.mu <- q.l(mu)(TruncQuantile)
- up0.mu <- q.l(mu)(TruncQuantile, lower.tail = FALSE)
- m0.mu <- median(mu); s0.mu <- IQR(mu)
- low1.mu <- m0.mu - IQR.fac * s0.mu
- up1.mu <- m0.mu + IQR.fac * s0.mu
- low.mu <- max(low0.mu,low1.mu); up.mu <- min(up0.mu,up1.mu)
-
-
- if(is(distr,"DiscreteDistribution"))
- x.seq <-support(distr)
- else
- {if(is(distr,"AbscontDistribution")){
- x.seq0 <- seq(low, up, length = N1)
- h0 <- diff(x.seq0[2:1])
- x.seq <- x.seq0[odd]
- }else{
- x.seq <- seq(low,up, length = N)
- }
- }
- if(is(mu,"DiscreteDistribution"))
- x.mu.seq <- support(mu)
- else
- {if(is(mu,"AbscontDistribution")){
- x.mu.seq0 <- seq(low.mu, up.mu, length = N1)
- h0.mu <- diff(x.mu.seq0[2:1])
- x.mu.seq <- x.mu.seq0[odd]
- }else{
- x.mu.seq <- seq(low.mu, up.mu, length = N)
- }
- }
-
- L2deriv <- L2deriv(L2Fam)[[1]]
-# y.seq <- sapply(x.seq, function(x) evalRandVar(L2deriv, x))
-# plot(x.seq[!is.na(y.seq)],y.seq ,type="l")
-
- ## are we working with a one-dim L2deriv or not?
-
- onedim <- (length(L2deriv at Map)==1)
-
-
- if(onedim){
- ## one-dim case
-
- ## Delta, formula (56), p. 133 [Ri:94]
- ## Ptheta- primitive function for Lambda
-
- if(is(distr,"AbscontDistribution")){
- Delta0x <- sapply(x.seq0, function(x)
- evalRandVar(L2deriv, x)) *
- d(distr)(x.seq0)
- Delta0 <- h0*.csimpsum(Delta0x)
- }else{
- L2x <- function(x,y) (x<=y)*evalRandVar(L2deriv, x)
- Delta0 <- sapply(x.seq, function(Y){ fct <- function(x) L2x(x,y=Y)
- return(E(object=distr, fun = fct))})
- }
- # print(Delta0)
- Delta1 <- approxfun(x.seq, Delta0, yleft = 0, yright = 0)
- if(is(distr,"DiscreteDistribution"))
- Delta <- function(x) Delta1(x) * (x %in% support(distr))
- else Delta <- function(x) Delta1(x)
- # print(Delta(x.seq))
- # print(Delta(rnorm(100)))
-
- ## J = Var_Ptheta Delta
- J1 <- E(object=distr, fun = Delta)
-# print(J1)
- Delta.0 <- function(x) Delta(x) - J1
- # print(Delta.0(x.seq))
- # print(Delta.0(r(distr)(100))^2)
- #J <- distrExIntegrate(function(x) d(distr)(x)*Delta.0(x)^2, lower=low, upper=up)
- J <- E(object=distr, fun = function(x) Delta.0(x)^2 )
-# print(J)
-
- ### CvM-IC phi
- phi <- function(x) Delta.0(x)/J
-
- ## integrand phi x Ptheta in formula (51) [ibid]
- phi1 <- function(x) phi(x) * p(distr)(x)
- psi1 <- E(object = mu, fun = phi1)
-
-
- ## obtaining IC psi (formula (51))
-
- if(is(mu,"AbscontDistribution")){
- phix <- function(x) phi(x)*d(mu)(x)
- psi0x <- sapply(rev(x.mu.seq0), phix)
- psi0 <- h0.mu*rev(.csimpsum(psi0x))
- }else{
- phixy <- function(x,y) (x<=y)*phi(y)
- psi0 <- sapply(x.mu.seq, function(X){ fct <- function(y) phixy(x=X,y=y)
- return(E(object=mu, fun = fct))})
- }
- # print(psi0)
- psi.1 <- approxfun(x.mu.seq, psi0, yleft = 0, yright = rev(psi0)[1])
- if(is(distr,"DiscreteDistribution"))
- psi <- function(x) (psi.1(x)-psi1) * (x %in% support(mu))
- else psi <- function(x) psi.1(x)-psi1
-
- E2 <- E(object=distr, fun = function(x) psi(x)^2)
- L2deriv <- L2Fam at L2deriv[[1]]
- ## E2 = Cov_mu (psi)
-
-# ### control: centering & standardization
- E1 <- E(object=distr, fun = psi )
- E3 <- E(object=distr, fun = function(x) psi(x)*evalRandVar(L2deriv, x))
- psi.0 <- function(x) psi(x) - E1
- psi.01 <- function(x) psi.0(x)/E3
- if(withplot)
- { dev.new() #windows()
- plot(x.seq, psi.01(x.seq),
- type = if(is(distr,"DiscreteDistribution")) "p" else "l")
- }
- E4 <- E(object=distr, fun = function(x) psi.01(x)^2)
- psi.01 <- EuclRandVariable(Map = list(psi.01), Domain = Reals())
-
-# print(list(E2,E4,E2-E4))
-
- }else{
-
- ## multivariate case
-
- Dim <- length(evalRandVar(L2deriv, 1))
-
- ## Delta, formula (56), p. 133 [Ri:94]
- ## Ptheta- primitive function for Lambda
-
- Map.Delta <- vector("list",Dim)
- # print("HLL")
- # print(x.seq0)
- for(i in 1:Dim)
- { if(is(distr,"AbscontDistribution")){
- #print(L2deriv at Map[[i]])
- fct0 <- sapply(x.seq0, L2deriv at Map[[i]]) *
- d(distr)(x.seq0)
- #print(fct0)
- Delta0 <- h0*.csimpsum(fct0)
- }else{
- fct0 <- function(x,y) L2deriv at Map[[i]](x)*(x<=y)
- Delta0 <- sapply(x.seq, function(Y){ fct <- function(x) fct0(x,y=Y)
- return(E(object=distr, fun = fct))})
- }
- #print(Delta0)
- Delta1 <- approxfun(x.seq, Delta0, yleft = 0, yright = 0)
- if(is(distr,"DiscreteDistribution"))
- Delta <- function(x) Delta1(x) * (x %in% support(distr))
- else Delta <- function(x) Delta1(x)
- Map.Delta[[i]] <- Delta
- env.i <- environment(Map.Delta[[i]]) <- new.env()
- assign("i", i, envir=env.i)
- assign("fct", fct, envir=env.i)
- assign("fct0", fct0, envir=env.i)
- assign("Delta", Delta, envir=env.i)
- assign("Delta0", Delta0, envir=env.i)
- assign("Delta1", Delta1, envir=env.i)
- if(withplot){
- dev.new()
- #windows()
- plot(x.seq, sapply(x.seq,Map.Delta[[i]]),
- type = if(is(distr,"DiscreteDistribution")) "p" else "l")
- }
-
- }
- Delta <- EuclRandVariable(Map = Map.Delta, Domain = Reals())
-
-
-
- ## J = Var_Ptheta Delta
- J1 <- E(object=distr, fun = Delta)
- Delta.0 <- Delta - J1
- J <- E(object=distr, fun = Delta.0 %*%t(Delta.0))
- ### CvM-IC phi
- phi <- as(solve(J)%*%Delta.0,"EuclRandVariable")
-
- ## integrand phi x Ptheta in formula (51) [ibid]
-
- Map.phi1 <- vector("list",Dim)
- for(i in 1:Dim)
- { Map.phi1[[i]] <- function(x) evalRandVar(phi,x)[i] * p(distr)(x)
- env.i <- environment(Map.phi1[[i]]) <- new.env()
- assign("i", i, envir=env.i)
- }
-
- phi1 <- EuclRandVariable(Map = Map.phi1, Domain = Reals())
- psi1 <- E(object=mu, fun = phi1)
-
-
- ## obtaining IC psi (formula (51))
- Map.psi <- vector("list",Dim)
- for(i in 1:Dim)
- { if(is(mu,"AbscontDistribution")){
- fct01 <- function(x) phi at Map[[i]](x)*d(mu)(x)
- fct0 <- sapply(rev(x.mu.seq0),fct01)
- phi0 <- h0.mu*rev(.csimpsum(fct0))
- }else{
- fct01 <- NULL
- fct0 <- function(x,y) evalRandVar(phi, y)[i]*(x<=y)
- phi0 <- sapply(x.mu.seq,
- function(X){
- fct <- function(y) fct0(x = X, y)
- return(E(object = mu, fun = fct))
- })
- }
-
- phi0a <- approxfun(x.mu.seq, phi0, yleft = 0, yright = rev(phi0)[1])
- env.i <- environment(phi1) <- new.env()
- assign("i", i, envir=env.i)
- if(is(distr,"DiscreteDistribution"))
- psi0 <- function(x) phi0a(x) * (x %in% support(mu))
- else psi0 <- function(x) phi0a(x)
-
- Map.psi[[i]] <- psi0
- env.i <- environment(Map.psi[[i]]) <- new.env()
- assign("i", i, envir=env.i)
- assign("fct", fct, envir=env.i)
- assign("fct0", fct0, envir=env.i)
- assign("psi0", psi0, envir=env.i)
- assign("phi0a", phi0a, envir=env.i)
- assign("phi0", phi0, envir=env.i)
- }
- psi <- EuclRandVariable(Map = Map.psi, Domain = Reals())
-
- E2 <- E(object=distr, fun = psi %*%t(psi))
- ## E2 = Cov_mu (psi)
-
- ### control: centering & standardization
- L2deriv <- L2Fam at L2deriv[[1]]
- E1 <- E(object=distr, fun = psi )
- E3 <- E(object=distr, fun = psi %*%t(L2deriv))
- psi.0 <- psi - E1
- psi.01 <- as(solve(E3)%*%psi.0,"EuclRandVariable")
- if(withplot)
- { for(i in 1:Dim)
- { dev.new()
- plot(x.mu.seq, sapply(x.mu.seq,psi.01 at Map[[i]]),
- type = if(is(distr,"DiscreteDistribution")) "p" else "l")
- }}
- E4 <- E(object=distr, fun = psi.01 %*%t(psi.01))
- }
- E4 <- PosSemDefSymmMatrix(E4)
-
- psi <- EuclRandVarList(psi.01)
- nms <- names(c(main(param(L2Fam)),nuisance(param(L2Fam))))
- dimnames(E4) = list(nms,nms)
- if(withpreIC) return(list(preIC=psi, Var=E4))
- else return(E4)
-}
-
-### examples:
-if(FALSE){
-P0 <- PoisFamily();.CvMMDCovariance(P0,par=ParamFamParameter("lambda",1), withplot=TRUE)
-B0 <- BinomFamily(size=8, prob=0.3);.CvMMDCovariance(B0,par=ParamFamParameter("",.3), withplot=TRUE)
-N0 <- NormLocationFamily();.CvMMDCovariance(N0,par=ParamFamParameter("",0), withplot=TRUE, N = 200)
-C0 <- L2LocationFamily(central=Cauchy());.CvMMDCovariance(C0,par=ParamFamParameter("",0), withplot=TRUE, N = 200)
-N1 <- NormScaleFamily(); re=.CvMMDCovariance(N1,par=ParamFamParameter("",1), withICwithplot=TRUE, N = 200)
-NS <- NormLocationScaleFamily();paramP <- ParamFamParameter(name = "locscale", main = c("loc"=0,"scale"=1),trafo = diag(2));
- .CvMMDCovariance(NS,par=paramP, withplot=TRUE, N = 100)
-cls <- CauchyLocationScaleFamily();.CvMMDCovariance(cls,par=ParamFamParameter("",0:1), withplot=TRUE, N = 200)
-Els <- L2LocationScaleFamily(loc = 0, scale = 1,
- name = "Laplace Location and scale family",
- centraldistribution = DExp(),
- LogDeriv = function(x) sign(x),
- FisherInfo = diag(2),
- trafo = diag(2))
-.CvMMDCovariance(Els,par=ParamFamParameter("",0:1), withplot=TRUE, N = 100)
-
-system.time(print(.CvMMDCovariance(P0,par=ParamFamParameter("lambda",1))))
-system.time(print(.CvMMDCovariance(B0,par=ParamFamParameter("",.3))))
-system.time(print(.CvMMDCovariance(N0,par=ParamFamParameter("",0), N = 100)))
-system.time(print(.CvMMDCovariance(C0,par=ParamFamParameter("",0), N = 100)))
-system.time(print(.CvMMDCovariance(N1,par=ParamFamParameter("",1), N = 100)))
-system.time(print(.CvMMDCovariance(NS,par=paramP, N = 100)))
-system.time(print(.CvMMDCovariance(cls,par=ParamFamParameter("",0:1), N = 100)))
-system.time(print(.CvMMDCovariance(Els,par=ParamFamParameter("",0:1), N = 100)))
-
-}
-
#------------------------------------
#### utilities copied from package distr v.2.6 svn-rev 943
#------------------------------------
Modified: branches/distr-2.8/pkg/distrMod/R/MCEstimator.R
===================================================================
--- branches/distr-2.8/pkg/distrMod/R/MCEstimator.R 2018-08-05 15:52:18 UTC (rev 1242)
+++ branches/distr-2.8/pkg/distrMod/R/MCEstimator.R 2018-08-05 15:55:36 UTC (rev 1243)
@@ -53,6 +53,7 @@
if(!is.null(asv)) argList <- c(argList, asvar.fct = asv)
if(!is.null(dots)) argList <- c(argList, dots)
+ argList <- c(argList, x = x)
## digesting the results of mceCalc
res <- do.call(.process.meCalcRes, argList)
Modified: branches/distr-2.8/pkg/distrMod/R/MDEstimator.R
===================================================================
--- branches/distr-2.8/pkg/distrMod/R/MDEstimator.R 2018-08-05 15:52:18 UTC (rev 1242)
+++ branches/distr-2.8/pkg/distrMod/R/MDEstimator.R 2018-08-05 15:55:36 UTC (rev 1243)
@@ -52,6 +52,7 @@
if(!is.null(dots)) argList <- c(argList, dots)
if(!validity.check %in% names(argList))
argList$validity.check <- TRUE
+ argList <- c(argList, x = x)
## digesting the results of mceCalc
res <- do.call(.process.meCalcRes, argList)
@@ -60,16 +61,32 @@
return(.checkEstClassForParamFamily(ParamFamily,res))
}
-CvMMDEstimator <- function(x, ParamFamily, paramDepDist = FALSE,
+CvMMDEstimator <- function(x, ParamFamily, muDatOrMod = c("Dat","Mod"),
+ paramDepDist = FALSE,
startPar = NULL, Infos,
trafo = NULL, penalty = 1e20,
validity.check = TRUE, asvar.fct = .CvMMDCovariance,
na.rm = TRUE, ..., .withEvalAsVar = TRUE){
- MDEstimator(x = x, ParamFamily = ParamFamily, distance = CvMDist,
+
+ muDatOrMod <- match.arg(muDatOrMod)
+ if(muDatOrMod=="Dat") {
+ distance0 <- CvMDist
+ estnsffx <- "(mu = emp. cdf)"
+ if(missing(asvar.fct)) asvar.fct <- .CvMMDCovarianceWithMux
+ }else{
+ distance0 <- CvMDist2
+ estnsffx <- "(mu = model distr.)"
+ if(missing(asvar.fct)) asvar.fct <- .CvMMDCovariance
+ }
+
+ res <- MDEstimator(x = x, ParamFamily = ParamFamily, distance = distance0,
paramDepDist = paramDepDist, startPar = startPar, Infos = Infos,
trafo = trafo, penalty = penalty, validity.check = validity.check,
asvar.fct = asvar.fct, na.rm = na.rm,
..., .withEvalAsVar = .withEvalAsVar)
+ res at name <- paste("Minimum CvM distance estimate", estnsffx)
+ res at estimate.call <- match.call()
+ return(res)
}
KolmogorovMDEstimator <- function(x, ParamFamily, paramDepDist = FALSE,
@@ -77,11 +94,13 @@
trafo = NULL, penalty = 1e20,
validity.check = TRUE, asvar.fct, na.rm = TRUE, ...,
.withEvalAsVar = TRUE){
- MDEstimator(x = x, ParamFamily = ParamFamily, distance = KolmogorovDist,
+ res <- MDEstimator(x = x, ParamFamily = ParamFamily, distance = KolmogorovDist,
paramDepDist = paramDepDist, startPar = startPar, Infos = Infos,
trafo = trafo, penalty = penalty, validity.check = validity.check,
asvar.fct = asvar.fct, na.rm = na.rm,
..., .withEvalAsVar = .withEvalAsVar)
+ res at estimate.call <- match.call()
+ return(res)
}
TotalVarMDEstimator <- function(x, ParamFamily, paramDepDist = FALSE,
@@ -89,11 +108,13 @@
trafo = NULL, penalty = 1e20,
validity.check = TRUE, asvar.fct, na.rm = TRUE, ...,
.withEvalAsVar = TRUE){
- MDEstimator(x = x, ParamFamily = ParamFamily, distance = TotalVarDist,
+ res <- MDEstimator(x = x, ParamFamily = ParamFamily, distance = TotalVarDist,
paramDepDist = paramDepDist, startPar = startPar, Infos = Infos,
trafo = trafo, penalty = penalty, validity.check = validity.check,
asvar.fct = asvar.fct, na.rm = na.rm,
..., .withEvalAsVar = .withEvalAsVar)
+ res at estimate.call <- match.call()
+ return(res)
}
HellingerMDEstimator <- function(x, ParamFamily, paramDepDist = FALSE,
@@ -101,10 +122,12 @@
trafo = NULL, penalty = 1e20,
validity.check = TRUE, asvar.fct, na.rm = TRUE, ...,
.withEvalAsVar = TRUE){
- MDEstimator(x = x, ParamFamily = ParamFamily, distance = HellingerDist,
+ res <- MDEstimator(x = x, ParamFamily = ParamFamily, distance = HellingerDist,
paramDepDist = paramDepDist, startPar = startPar, Infos = Infos,
trafo = trafo, penalty = penalty, validity.check = validity.check,
asvar.fct = asvar.fct, na.rm = na.rm,
..., .withEvalAsVar = .withEvalAsVar)
+ res at estimate.call <- match.call()
+ return(res)
}
Modified: branches/distr-2.8/pkg/distrMod/R/MLEstimator.R
===================================================================
--- branches/distr-2.8/pkg/distrMod/R/MLEstimator.R 2018-08-05 15:52:18 UTC (rev 1242)
+++ branches/distr-2.8/pkg/distrMod/R/MLEstimator.R 2018-08-05 15:55:36 UTC (rev 1243)
@@ -48,6 +48,7 @@
if(!is.null(asv)) argList <- c(argList, asvar.fct = asv)
if(!is.null(dots)) argList <- c(argList, dots)
+ argList <- c(argList, x = x)
## digesting the results of mceCalc
res <- do.call(what = ".process.meCalcRes", args = argList)
Modified: branches/distr-2.8/pkg/distrMod/R/SimpleL2ParamFamilies.R
===================================================================
--- branches/distr-2.8/pkg/distrMod/R/SimpleL2ParamFamilies.R 2018-08-05 15:52:18 UTC (rev 1242)
+++ branches/distr-2.8/pkg/distrMod/R/SimpleL2ParamFamilies.R 2018-08-05 15:55:36 UTC (rev 1243)
@@ -30,7 +30,7 @@
prob.0 <- main(param)
fct <- function(x){}
body(fct) <- substitute({ (x-size*prob.1)/(prob.1*(1-prob.1)) },
- list(size = size, prob.0 = prob.1))
+ list(size = size, prob.1 = prob.0))
return(fct)}
L2derivSymm <- FunSymmList(OddSymmetric(SymmCenter = size*prob))
L2derivDistr <- UnivarDistrList((distribution - size*prob)/(prob*(1-prob)))
@@ -45,7 +45,7 @@
FisherInfo <- FisherInfo.fct(param.0)
res <- L2ParamFamily(name = name, distribution = distribution,
- distrSymm = distrSymm, param = param, modifyParam = modifyParam,
+ distrSymm = distrSymm, param = param.0, modifyParam = modifyParam,
props = props, L2deriv.fct = L2deriv.fct, L2derivSymm = L2derivSymm,
L2derivDistr = L2derivDistr, L2derivDistrSymm = L2derivDistrSymm,
FisherInfo.fct = FisherInfo.fct, FisherInfo = FisherInfo,
@@ -97,7 +97,7 @@
FisherInfo <- FisherInfo.fct(param.0)
res <- L2ParamFamily(name = name, distribution = distribution,
- distrSymm = distrSymm, param = param, modifyParam = modifyParam,
+ distrSymm = distrSymm, param = param.0, modifyParam = modifyParam,
props = props, L2deriv.fct = L2deriv.fct, L2derivSymm = L2derivSymm,
L2derivDistr = L2derivDistr, L2derivDistrSymm = L2derivDistrSymm,
FisherInfo.fct = FisherInfo.fct, FisherInfo = FisherInfo,
@@ -155,7 +155,7 @@
FisherInfo <- FisherInfo.fct(param.0)
res <- L2ParamFamily(name = name, distribution = distribution,
- distrSymm = distrSymm, param = param, modifyParam = modifyParam,
+ distrSymm = distrSymm, param = param.0, modifyParam = modifyParam,
props = props, L2deriv.fct = L2deriv.fct, L2derivSymm = L2derivSymm,
L2derivDistr = L2derivDistr, L2derivDistrSymm = L2derivDistrSymm,
FisherInfo.fct = FisherInfo.fct, FisherInfo = FisherInfo,
@@ -229,7 +229,7 @@
FisherInfo <- FisherInfo.fct(param.0)
res <- L2ParamFamily(name = name, distribution = distribution,
- distrSymm = distrSymm, param = param, modifyParam = modifyParam,
+ distrSymm = distrSymm, param = param.0, modifyParam = modifyParam,
props = props, L2deriv.fct = L2deriv.fct, L2derivSymm = L2derivSymm,
L2derivDistr = L2derivDistr, L2derivDistrSymm = L2derivDistrSymm,
FisherInfo.fct = FisherInfo.fct, FisherInfo = FisherInfo,
@@ -318,7 +318,7 @@
FisherInfo <- FisherInfo.fct(param.0)
res <- L2ParamFamily(name = name, distribution = distribution,
- distrSymm = distrSymm, param = param, modifyParam = modifyParam,
+ distrSymm = distrSymm, param = param.0, modifyParam = modifyParam,
props = props, L2deriv.fct = L2deriv.fct, L2derivSymm = L2derivSymm,
L2derivDistr = L2derivDistr, L2derivDistrSymm = L2derivDistrSymm,
FisherInfo.fct = FisherInfo.fct, FisherInfo = FisherInfo,
@@ -394,7 +394,7 @@
L2Fam at name <- name
L2Fam at distribution <- distribution
L2Fam at distrSymm <- distrSymm
- L2Fam at param <- param
+ L2Fam at param <- param.0
L2Fam at modifyParam <- modifyParam
L2Fam at props <- props
L2Fam at L2deriv.fct <- L2deriv.fct
@@ -406,7 +406,7 @@
L2Fam at startPar <- startPar
L2Fam at makeOKPar <- makeOKPar
L2Fam at scaleshapename <- c("scale"="scale","shape"="shape")
-
+
L2deriv <- EuclRandVarList(RealRandVariable(L2deriv.fct(param.0),
Domain = Reals()))
@@ -427,7 +427,6 @@
return(L2Fam)
}
-(G1 <- GammaFamily())
##################################################################
## Beta family :: new 08/08 P.R.
@@ -483,7 +482,7 @@
FisherInfo <- FisherInfo.fct(param.0)
res <- L2ParamFamily(name = name, distribution = distribution,
- distrSymm = distrSymm, param = param, modifyParam = modifyParam,
+ distrSymm = distrSymm, param = param.0, modifyParam = modifyParam,
props = props, L2deriv.fct = L2deriv.fct, L2derivSymm = L2derivSymm,
L2derivDistr = L2derivDistr, L2derivDistrSymm = L2derivDistrSymm,
FisherInfo.fct = FisherInfo.fct, FisherInfo = FisherInfo,
Added: branches/distr-2.8/pkg/distrMod/R/asCvMVarianceQtl.R
===================================================================
--- branches/distr-2.8/pkg/distrMod/R/asCvMVarianceQtl.R (rev 0)
+++ branches/distr-2.8/pkg/distrMod/R/asCvMVarianceQtl.R 2018-08-05 15:55:36 UTC (rev 1243)
@@ -0,0 +1,780 @@
+.CvMMDCovarianceWithMux <- function(L2Fam, param, withplot = FALSE, withpreIC = FALSE,
+ N = 400, rel.tol=.Machine$double.eps^0.3,
+ TruncQuantile = getdistrOption("TruncQuantile"),
+ IQR.fac = 15, ..., x=NULL){
+ mu <- distribution(L2Fam)
+ if(!is.null(x)) mu <- DiscreteDistribution(x)
+ .CvMMDCovariance(L2Fam=L2Fam, param=param, mu=mu,
+ withplot = withplot, withpreIC = withpreIC,
+ N = N, rel.tol=rel.tol, TruncQuantile = TruncQuantile,
+ IQR.fac = IQR.fac, ...)
+}
+
+CvMDist2 <- function(e1,e2,... ) CvMDist(e1, e2, mu = e2, ...)
+
+### 20180805: new function to compute asCov of CvM-MDE
+# which for the primitive functions uses integration on [0,1]
+# via quantile transformation
+.CvMMDCovariance<- function(L2Fam, param, mu = distribution(L2Fam),
+ withplot = FALSE, withpreIC = FALSE,
+ N = 400, rel.tol=.Machine$double.eps^0.3,
+ TruncQuantile = getdistrOption("TruncQuantile"),
+ IQR.fac = 15,
+ ...){
+ # preparations:
+
+ dotsInt <- list(...)
+ dotsInt[["f"]] <- NULL
+ dotsInt[["lower"]] <- NULL
+ dotsInt[["upper"]] <- NULL
+ dotsInt[["stop.on.error"]] <- NULL
+ dotsInt[["distr"]] <- NULL
+
+ N.1 <- round(0.2*N)
+ N.3 <- N.1
+ N.2 <- N-N.1-N.3
+
+ N1 <- 2*N+1
+ N1.1 <- 2*N.1+1
+ N1.2 <- 2*N.2+1
+ N1.3 <- 2*N.3+1
+ odd <- (1:N1)%%2==1
+ odd.1 <- (1:N1.1)%%2==1
+ odd.2 <- (1:N1.2)%%2==1
+ odd.3 <- (1:N1.3)%%2==1
+
+ param0 <- L2Fam at param
+ dim0 <- dimension(param0)
+
+ paramP <- param0
+ paramP at main <- main(param)
+ paramP at trafo <- diag(dim0)
+ L2Fam <- modifyModel(L2Fam, paramP)
+
+ distr <- L2Fam at distribution
+
+ ### get a sensible integration range:
+ low <- TruncQuantile
+ up <- 1-TruncQuantile
+
+
+ if(is(distr,"DiscreteDistribution"))
+ x.seq <-support(distr)
+ else
+ {if(is(distr,"AbscontDistribution")){
+
+ ## split up the integration range into
+ # .1 = lower tail, .2 mid range, .3 upper tail
+
+ x.seq0 <- seq(0, 1, length = N1)
+ h0 <- diff(x.seq0[1:2])
+ x.seq0.1 <- seq(0, 1, length = N1.1)
+ h0.1 <- diff(x.seq0.1[1:2])
+ x.seq0.2 <- seq(0, 1, length = N1.2)
+ h0.2 <- diff(x.seq0.2[1:2])
+ x.seq <- x.seq0[odd]
+ x.seq.1 <- low+(1-low)*x.seq0.1/100
+ x.seq.3 <- 1-rev(x.seq.1)
+ x.seq.la <- rev(x.seq.1)[1]
+ del <- 1-2*(x.seq.la+(h0.1/100+h0.2)/2)
+ x.seq.2l <- x.seq.la+(h0.1/100+h0.2)/2+del*x.seq0.2
+ x.seq.2r <- 1-rev(x.seq.2l)
+ x.seq.2 <- (x.seq.2l+x.seq.2r)/2
+ x.seq.a <- c(x.seq.1[odd.1],x.seq.2[odd.2],x.seq.3[odd.3])
+# x.seq.b <- c(x.seq.1,x.seq.2,x.seq.3)
+# iN.1 <- 1:N1.1
+# iN.2 <- N1.1+(1:N1.2)
+# iN.3 <- N1.1+N1.2+(1:N1.3)
+# riN.3 <- 1:N1.3
+# riN.2 <- N1.3+1:N1.2
+# riN.1 <- N1.3+N1.2+1:N1.1
+ }else{
+ x.seq <- seq(low,up, length = N)
+ }
+ }
+ if(is(mu,"DiscreteDistribution"))
+ x.mu.seq <- support(mu)
+ else
+ {if(is(mu,"AbscontDistribution")){
+ x.mu.seq0 <- x.seq0
+ h0.mu <- h0
+ x.mu.seq <- x.seq
+ x.mu.seq.1 <- x.seq.1
+ x.mu.seq.2 <- x.seq.2
+ x.mu.seq.3 <- x.seq.3
+ x.mu.seq.a <- x.seq.a
+# x.mu.seq.b <- x.seq.b
+# iN.mu.1 <- iN.1
+# iN.mu.2 <- iN.2
+# iN.mu.3 <- iN.3
+# riN.mu.1 <- riN.1
+# riN.mu.2 <- riN.2
+# riN.mu.3 <- riN.3
+ }else{
+ x.mu.seq <- seq(low, up, length = N)
+ }
+ }
+
+ L2deriv.0 <- L2deriv(L2Fam)[[1]]
+# y.seq <- sapply(x.seq, function(x) evalRandVar(L2deriv, x))
+# plot(x.seq[!is.na(y.seq)],y.seq ,type="l")
+
+ ## are we working with a one-dim L2deriv or not?
+
+ onedim <- (length(L2deriv.0 at Map)==1)
+
+
+ myint <- function(f,...){
+ distrExIntegrate(f=f, lower=0, upper=1,
+ stop.on.error=FALSE, distr=Unif(), ...)
+ }
+
+ if(onedim){
+ ## one-dim case
+
+ ## Delta, formula (56), p. 133 [Ri:94]
+ ## Ptheta- primitive function for Lambda
+
+ if(is(distr,"AbscontDistribution")){
+ fqx <- function(x){qx <- q.l(distr)(x)
+ return(sapply(qx,function(y)evalRandVar(L2deriv.0, y)))
+ }
+ #Delta0x <- sapply(x.seq.b,fqx)
+ #Delta0x.1 <- Delta0x[iN.1]
+ #Delta0x.2 <- Delta0x[iN.2]
+ #Delta0x.3 <- Delta0x[iN.3]
+ Delta0x.1 <- sapply(x.seq.1,fqx)
+ Delta0x.2 <- sapply(x.seq.2,fqx)
+ Delta0x.3 <- sapply(x.seq.3,fqx)
+ Delta0.1 <- h0/100*.csimpsum(Delta0x.1)
+ Delta0.2 <- rev(Delta0.1)[1]+h0*.csimpsum(Delta0x.2)
+ Delta0.3 <- rev(Delta0.2)[1]+h0/100*.csimpsum(Delta0x.3)
+ Delta0 <- c(Delta0.1,Delta0.2,Delta0.3)
+ Delta1.q <- approxfun(x.seq.a, Delta0, yleft = 0, yright = 0)
+ J1 <- do.call(myint, c(list(f=Delta1.q), dotsInt))
+ Delta.0 <- function(x) Delta1.q(p(distr)(x))-J1
+ J <- do.call(myint, c(list(f=function(x) (Delta1.q(x)-J1)^2),dotsInt))
+ }else{
+ L2x <- function(x,y) (x<=y)*evalRandVar(L2deriv.0, x)
+ Delta0 <- sapply(x.seq, function(Y){ fct <- function(x) L2x(x,y=Y)
+ return(E(object=distr, fun = fct))})
+ Delta1 <- approxfun(x.seq, Delta0, yleft = 0, yright = 0)
+ Delta <- Delta1
+ if(is(distr,"DiscreteDistribution"))
+ Delta <- function(x) Delta1(x) * (x %in% support(distr))
+ J1 <- E(object=distr, fun = Delta)
+ Delta.0 <- function(x) Delta(x) - J1
+ J <- E(object=distr, fun = function(x) Delta.0(x)^2 )
+ }
+
+ ### CvM-IC phi
+ phi <- function(x) Delta.0(x)/J
+
+ ## obtaining IC psi (formula (51))
+
+ if(is(mu,"AbscontDistribution")){
+ ## integrand phi x Ptheta in formula (51) [ibid]
+ phi1.q <- function(s){qs <- q.l(mu)(s)
+ return(phi(qs)*p(distr)(qs)) }
+ psi1 <- do.call(myint, c(list(f=phi1.q),dotsInt))
+
+ phiqx <- function(x){qx <- q.l(mu)(x)
+ return(phi(qx))}
+ #psi0qx <- sapply(rev(x.mu.seq.b), phiqx)
+ #psi0qx.1 <- psi0qx[riN.mu.1]
+ #psi0qx.2 <- psi0qx[riN.mu.2]
+ #psi0qx.3 <- psi0qx[riN.mu.3]
+ psi0qx.1 <- sapply(rev(x.mu.seq.1), phiqx)
+ psi0qx.2 <- sapply(rev(x.mu.seq.2), phiqx)
+ psi0qx.3 <- sapply(rev(x.mu.seq.3), phiqx)
+ psi0q.3 <- h0.mu/100*rev(.csimpsum(psi0qx.3))
+ psi0q.2 <- psi0q.3[1]+h0.mu*rev(.csimpsum(psi0qx.2))
+ psi0q.1 <- psi0q.2[1]+h0.mu/100*rev(.csimpsum(psi0qx.1))
+
+ psi0q <- c(psi0q.1,psi0q.2,psi0q.3)
+ psi.q1 <- approxfun(x.mu.seq.a, psi0q, yleft = 0, yright = rev(psi0q)[1])
+ psi <- function(x) psi.q1(p(mu)(x))-psi1
+ }else{
+ ## integrand phi x Ptheta in formula (51) [ibid]
+ phi1 <- function(x) phi(x) * p(distr)(x)
+ psi1 <- E(object = mu, fun = phi1)
+
+ phixy <- function(x,y) (x<=y)*phi(y)
+ psi0 <- sapply(x.mu.seq, function(X){ fct <- function(y) phixy(x=X,y=y)
+ return(E(object=mu, fun = fct))})
+ psi.1 <- approxfun(x.mu.seq, psi0, yleft = 0, yright = rev(psi0)[1])
+ psi <- function(x) psi.1(x)-psi1
+ if(is(distr,"DiscreteDistribution"))
+ psi <- function(x) (psi.1(x)-psi1) * (x %in% support(mu))
+ }
+ # print(psi0)
+ if(is(distr,"AbscontDistribution")){
+ psi.q <- function(x){qx <- q.l(distr)(x); return(psi(qx))}
+ E2 <- do.call(myint, c(list(f=function(x)psi.q(x)^2),dotsInt))
+ E1 <- do.call(myint, c(list(f=psi.q),dotsInt))
+ E3 <- do.call(myint, c(list(f=function(x){
+ qx <- q.l(distr)(x)
+ L2qx <- sapply(qx,function(y)
+ evalRandVar(L2deriv.0, y))
+ return(psi(qx)*L2qx)
+ }), dotsInt))
+ psi.01 <- function(x) (psi(x)-E1)/E3
+ E4 <- do.call(myint, c(list(f=function(x) (psi.q(x)-E1)^2/E3^2),dotsInt))
+ }else{
+ E2 <- E(object=distr, fun = function(x) psi(x)^2)
+ L2x <- function(x,y) (x<=y)*evalRandVar(L2deriv.0, x)
+ E1 <- E(object=distr, fun = psi )
+ E3 <- E(object=distr, fun = function(x) psi(x)*evalRandVar(L2deriv.0, x))
+ psi.0 <- function(x) psi(x) - E1
+ psi.01 <- function(x) psi.0(x)/E3
+ E4 <- E(object=distr, fun = function(x) psi.01(x)^2)
+ }
+ ## E2 = Cov_mu (psi)
+
+# ### control: centering & standardization
+ if(withplot)
+ { dev.new() #windows()
+ x0.seq <- x.seq
+ if(is(distr,"AbscontDistribution")) x0.seq <- q.l(distr)(x.seq)
+ plot(x0.seq, psi.01(x0.seq),
+ type = if(is(distr,"DiscreteDistribution")) "p" else "l")
+ }
+ psi.01 <- EuclRandVariable(Map = list(psi.01), Domain = Reals())
+
+
+ }else{
+
+ ## multivariate case
+
+ Dim <- length(evalRandVar(L2deriv.0, 1))
+
+ ## Delta, formula (56), p. 133 [Ri:94]
+ ## Ptheta- primitive function for Lambda
+
+ Map.Delta <- vector("list",Dim)
+
+ for(i in 1:Dim)
+ { if(is(distr,"AbscontDistribution")){
+ #fct0.q <- sapply(x.seq.b, function(x){qx <- q.l(distr)(x); return(L2deriv.0 at Map[[i]](qx))})
+ #fct0.q1 <- fct0.q[iN.1]
+ #fct0.q2 <- fct0.q[iN.2]
+ #fct0.q3 <- fct0.q[iN.3]
+ fct0.q1 <- sapply(x.seq.1, function(x){qx <- q.l(distr)(x); return(L2deriv.0 at Map[[i]](qx))})
+ fct0.q2 <- sapply(x.seq.2, function(x){qx <- q.l(distr)(x); return(L2deriv.0 at Map[[i]](qx))})
+ fct0.q3 <- sapply(x.seq.3, function(x){qx <- q.l(distr)(x); return(L2deriv.0 at Map[[i]](qx))})
+ #print(fct0)
+ Delta0.q1 <- h0/100*.csimpsum(fct0.q1)
+ Delta0.q2 <- rev(Delta0.q1)[1]+h0*.csimpsum(fct0.q2)
+ Delta0.q3 <- rev(Delta0.q2)[1]+h0/100*.csimpsum(fct0.q3)
+ Delta0.q <- c(Delta0.q1,Delta0.q2,Delta0.q3)
+ Delta1.q <- approxfun(x.seq.a, Delta0.q, yleft = 0, yright = 0)
+ Delta <- function(x) Delta1.q(p(distr)(x))
+ Map.Delta[[i]] <- Delta
+ env.i <- environment(Map.Delta[[i]]) <- new.env()
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/distr -r 1243
More information about the Distr-commits
mailing list