wthom.test <- function(fx, tf, x, wt, theta){

  #
  # description:
  # Computes the weighted homogeneity test statistic for a test of the
  # number of components in a mixture.
  #
  # arguments: 
  # fx - a function that computes the component density for the mixture.
  # tf -  a function that computes the homogeneity test function t(x,theta).
  # x - the vector or the matrix of observations. For a matrix, individual
  #     observation vectors should be stored in columns.
  # wt - the vector of weight parameters.
  # theta - a vector or matrix of component parameters. The jth column should
  #         give the parameters for the jth component of the mixture.
  #
  # details:
  # The function fx should be of the form fx(d,x,theta) and return the vector
  # of component densities (d = 0) for the vector (matrix) x of all observations
  # and single parameter theta. If d=1, it should return the matrix of partial
  # derivatives. Each row should give the vector of partial derivatives for the
  # corresponding x value.
  #
  # The function tx should be of the form tx(x,theta) and return the
  # vector of homogeneity test functions t(x,theta) for the vector (matrix)
  # x of all observations and a single parameter theta.
  #
  # value:
  # tstat - the test statistics for each of the components.
  # se - the corresponding standard error estimates.
  #
  # examples:
  # mu <- c(-2.752136, -0.416129)
  # v <- c(0.540891, 3.361862)
  # wt <- c(0.198672, 0.801328)
  # theta <- rbind(mu, v)
  # ts <- wthom.test(wthom.norm.fx, wthom.norm.tf, x, wt, theta)
  # pvalue <- 2*pnorm(-abs(ts$ts/ts$se))
  #

  # dimensionality parameters
  if (!is.matrix(theta)) theta <- matrix(theta, nrow = 1)
  m <- length(wt)
  p <- dim(theta)[1]
  if (m != dim(theta)[2])
    stop("second dimension of theta should equal length of wt")
  if (is.matrix(x))
    n <- dim(x)[2]
  else
    n <- length(x)
  di <- m-1+p*m

  tw <- NULL # matrix of density values 
  for(j in 1:m) tw <- cbind(tw, fx(0, x, theta[,j]))
    
  fm <- tw %*% wt # mixture density vector
  tw <- tw / c(fm)    # weight matrix
  
  u <- tw[,1:(m-1)] - tw[,m] # matrix of score vectors
  for(j in 1:m) u <- cbind(u, wt[j]*fx(1, x, theta[,j])/c(fm))

  afi <- t(u) %*% u / n  # approximate Fisher info
  
  # matrix of t(x,theta) w(x,theta)
  for(j in 1:m) tw[,j] <- tw[,j]*tf(x, theta[,j])

  tstat <- apply(tw, 2, mean) # un-normalized test statistics
  se <- apply(tw**2, 2, mean) # standard error

  corr <- NULL # correction for estimation to standard error
  for(j in 1:m){
    eutw <- apply(u*tw[,j], 2, mean)
    corr <- c(corr, t(eutw) %*% solve(afi, eutw))
  }

  se <- sqrt((se - corr)/n)
  return(tstat, se)
}

wthom.norm.fx <- function(d, x, theta){

  #
  # description:
  # computes the normal density (d=0) or partial derivative matrix.
  #
  # see also:
  # wthom.test
  
  mu <- theta[1]
  v <- theta[2]
  f <- exp( -(x - mu)**2/(2.0*v) )/sqrt(v) # note un-normalized
  if (d == 0) return(f)
  if (d == 1){
    fd <- f*(x - mu)/v
    fd <- cbind(fd, f * 0.5*(-1.0/v + ((x - mu)/v)**2))
    return(fd)
  }
}

wthom.norm.tf <- function(x, theta){

  # description:
  # the homogeneity test statistic t(x,theta) for use with normal distribution
  #      t(x,theta) = {|x-mu|/v}**3 - 4/root(2pi)
  #
  f <- abs(x-theta[1])/sqrt(theta[2])
  f <- f**3 - 4.0/sqrt(2.0*pi)
}

wthom.pois.fx <- function(d, x, theta){

  #
  # description:
  # computes the Poisson density (d=0) or derivative (d=1)_.
  #
  # see also:
  # wthom.test
  #

  if (d == 0) f <- exp(-theta)*theta**x #note: unormalized
  if (d == 1){
    if (theta > 0.0){
      f <- -exp(-theta)*theta**x
      f <- f + ifelse(x > 0, x*exp(-theta)*theta**(x-1), 0)
    }
    else{
      f <- ifelse(x > 1, 0.0, 1.0)
      f <- ifelse(x > 0, f, -1.0)
    }
  }

  return(f)
}

wthom.pois.tf <- function(x, theta){

  # description:
  # the homogeneity test statistic t(x,theta) for use with Poisson distribution
  #      t(x,theta) = (x - theta)**2 - x

  return((x-theta)**2 - x)
}


wthom.binom.fx <- function(d, x, theta){

  #
  # description:
  # computes the binomial density (d=0) or derivative (d=1).
  #
  # arguments:
  # x - the observations are assumed across columns. x[,i] gives the ith
  # observation: the number of successes, x[1,i], out of x[2,i] trials.
  #
  # see also:
  # wthom.test
  #

  ix <- x[1,]
  ni <- x[2,]
  if (d == 0)
    f <- theta**ix * (1-theta)**(ni - ix) #note: unnormalized
  if (d > 0){
    if (theta > 0.0 & theta < 1.0){
      f <- ix * theta**(ix - 1) * (1-theta)**(ni - ix) -
        (ni - ix) * theta**ix * (1-theta)**(ni - ix - 1)
    }
    if (theta >= 1.0){
      f <- ifelse(ix == ni, ni, -1.0)
      f <- ifelse(ix >= ni - 1, f, 0.0)
    }
    if (theta <= 0.0){
      f <- ifelse(ix == 0, ni, 1.0)
      f <- ifelse(ix >= 1, f, 0.0)
    }
  }

  return(f)
}

wthom.binom.tf <- function(x, theta){

  # description:
  # the homogeneity test statistic t(x,theta) for use with binomial distribution.
  #      
  # arguments:
  # x - the observations are assumed across columns. x[,i] gives the ith
  # observation: the number of successes, x[1,i], out of x[2,i] trials.
  
  ix <- x[1,]
  ni <- x[2,]
  iy <- ni - ix
  f <- (ix**2 - ix)/(theta) + (iy**2 - iy)/(1.0 - theta) - (ni**2 - ni)

  return(f)
}