.devel/R-old/p2.htest.approx.R

In Rexamine/agop: Aggregation Operators and Preordered Sets

## This file is part of the 'agop' library.
##
## Copyright 2013 Marek Gagolewski, Anna Cena
##
## Parts of the code are taken from the 'CITAN' R package by Marek Gagolewski
##
## 'agop' is free software: you can redistribute it and/or modify
## it under the terms of the GNU Lesser General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## 'agop' is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU Lesser General Public License for more details.
##
## You should have received a copy of the GNU Lesser General Public License
## along with 'agop'. If not, see <http://www.gnu.org/licenses/>.


invisible(NULL)

# #' Performs asymptotic (approximate) \eqn{h}-test for equality of shape parameters
# #' of two samples from the Pareto type-II distributions with known
# #' and equal scale parameters, \eqn{s>0}.
# #'
# #' Given two equal-sized samples \eqn{X=(X_1,...,X_n)} i.i.d. \eqn{P2(k_x,s)}
# #' and \eqn{Y=(Y_1,...,Y_m)} i.i.d. \eqn{P2(k_y,s)}
# #' this test verifies the null hypothesis
# #' \eqn{H_0: k_x=k_y}
# #' against two-sided or one-sided alternatives, depending
# #' on the value of \code{alternative}.
# #' It bases on a test statistic that is a function of {H(Y)-H(X)},
# #' where \eqn{H} denotes Hirsch's \eqn{h}-index (see \code{\link{index.h}}).
# #' This statistic approximately has asymptotically standardized normal distribution under \eqn{H_0}.
# #'
# #' Note that for \eqn{k_x < k_y}, then \eqn{X} dominates \eqn{Y} stochastically.
# #'
# #' @title Two-sample asymptotic h-test for equality of shape parameters for Type II-Pareto distributions with known common scale parameter
# #' @param x an n-element non-negative numeric vector of data values.
# #' @param y an n-element non-negative numeric vector of data values.
# #' @param s scale parameter, \eqn{s>0}.
# #' @param alternative indicates the alternative hypothesis and must be one of "two.sided" (default), "less", or "greater".
# #' @param significance significance level, \eqn{0<}\code{significance}\eqn{<1} or \code{NULL}. See Value for details.
# #' @return
# #' If \code{significance} is not \code{NULL}, then
# #' the list of class \code{power.htest} with the following components is passed as a result:
# #' \tabular{ll}{
# #' \code{statistic} \tab	the value of the test statistic.\cr
# #' \code{result} \tab	either FALSE (accept null hypothesis) or TRUE (reject).\cr
# #' \code{alternative} \tab	a character string describing the alternative hypothesis.\cr
# #' \code{method} \tab	a character string indicating what type of test was performed.\cr
# #' \code{data.name} \tab	a character string giving the name(s) of the data.\cr
# #' }
# #' Otherwise, the list of class \code{htest} with the following components is passed as a result:
# #' \tabular{ll}{
# #' \code{statistic} \tab	the value of the test statistic.\cr
# #' \code{p.value} \tab	the p-value of the test.\cr
# #' \code{alternative} \tab	a character string describing the alternative hypothesis.\cr
# #' \code{method} \tab	a character string indicating what type of test was performed.\cr
# #' \code{data.name} \tab	a character string giving the name(s) of the data.\cr
# #' }
# #' @export
# #' @seealso \code{\link{dpareto2}}, \code{\link{pareto2.goftest}}, \code{\link{pareto2.ftest}}, \code{\link{pareto2.htest}}, \code{\link{index.h}}
# #' @references
# #' Gagolewski M., Grzegorzewski P., S-Statistics and Their Basic Properties, In: Borgelt C. et al (Eds.),
# #' Combining Soft Computing and Statistical Methods in Data Analysis, Springer-Verlag, 2010, 281-288.\cr
# pareto2.htest.approx <- function(x, y, s, alternative = c("two.sided", "less", "greater"), significance=NULL)
# {
# 	alternative <- match.arg(alternative);
# 	DNAME <- deparse(substitute(x));
# 	DNAME <- paste(DNAME, "and", deparse(substitute(y)));
#
# 	if (mode(s) != "numeric" || length(s) != 1 || s <= 0) stop("'s' should be > 0");
#
# 	if (mode(x) != "numeric" || mode(y) != "numeric") stop("non-numeric data given");
#
# 	x <- x[!is.na(x)];
# 	n <- length(x);
# 	if (n < 1L || any(x<0)) stop("incorrect 'x' data");
#
# 	y <- y[!is.na(y)];
# 	if (length(y) < 1L || any(y<0)) stop("incorrect 'y' data");
#
# 	if (length(y) != n) stop("non-equal-sized vectors given on input");
#
#
#
# 	METHOD <- "Two-sample asymptotic h-test for equality of shape parameters for Type II-Pareto distributions with known common scale parameter";
#
# 	nm_alternative <- switch(alternative, two.sided = "two-sided",
# 			less = "kx < ky",
# 			greater = "kx > ky");
#
#
# 	# -----------------------------------------------------------------------
#
# 	HYn <- index.h(y,disable.check=TRUE)/n;
# 	HXn <- index.h(x,disable.check=TRUE)/n;
#
# 	kappa    <- function(x) { pmax(0,pmin(1,x))*n; }
#
# # 	if (HXn < 1e-9)
# # 	{
# # 		gprimex <- 0.0;
# # 	} else {
# 		kx <- uniroot(function(k,s,targetrho,kappa)
# 				{
# 					1-ppareto2(kappa(targetrho),k,s)-targetrho;
# 				}, c(1e-15,1e15), s, HXn, kappa, tol=1e-20)$root;
# 		gprimex <- dpareto2(HXn*n,kx,s)*n;
# # 	}
#
# # 	if (HYn < 1e-9)
# # 	{
# # 		gprimey <- 0.0;
# # 	} else {
# 		ky <- uniroot(function(k,s,targetrho,kappa)
# 				{
# 					1-ppareto2(kappa(targetrho),k,s)-targetrho;
# 				}, c(1e-15,1e15), s, HYn, kappa, tol=1e-20)$root;
# 		gprimey <- dpareto2(HYn*n,ky,s)*n;
# # 	}
#
# 	sigmax2 <- HXn*(1-HXn)/n/(1+gprimex)^2;
# 	sigmay2 <- HYn*(1-HYn)/n/(1+gprimey)^2;
#
#
# 	STATISTIC <- (HYn-HXn)/sqrt(sigmax2+sigmay2);
#
# 	names(STATISTIC) <- "T";
#
# 	# -----------------------------------------------------------------------
#
# 	if (!is.null(significance))
# 	{
# 		if (length(significance) != 1 || significance <= 0 || significance >=1)
# 			stop("incorrect 'significance'");
#
# 		if (significance > 0.2) warning("'significance' is possibly incorrect");
#
# 		RESULT <- ifelse(alternative == "two.sided", (STATISTIC<qnorm(significance*0.5) || STATISTIC>qnorm(1-significance*0.5)),
# 		          ifelse(alternative == "greater",    STATISTIC>qnorm(1-significance),
# 		                                              STATISTIC<qnorm(significance)));
#
# 		RVAL <- list(statistic = STATISTIC, result = RESULT, alternative = nm_alternative,
# 			method = METHOD, data.name = DNAME);
# 		class(RVAL) <- "power.htest";
# 		return(RVAL);
# 	} else {
#
# 		PVAL <- ifelse(alternative == "two.sided", (0.5-abs(pnorm(STATISTIC)-0.5))*2,
# 		        ifelse(alternative == "greater",   1-pnorm(STATISTIC),
# 		                                             pnorm(STATISTIC)));
#
#
# 		RVAL <- list(statistic = STATISTIC, p.value = PVAL, alternative = nm_alternative,
# 			method = METHOD, data.name = DNAME);
# 		class(RVAL) <- "htest";
# 		return(RVAL);
# 	}
# }