emplik2 (version 1.32)

el2.test.wtm: Computes maximum-likelihood probability jumps for multiple mean-type hypotheses, based on two independent uncensored samples

Description

This function computes the maximum-likelihood probability jumps for multiple mean-type hypotheses, based on two samples that are independent, uncensored, and weighted. The target function for the maximization is the constrained log(empirical likelihood) which can be expressed as, $$ \sum_{dx_i=1} wx_i \log{\mu_i} + \sum_{dy_j=1} wy_j \log{\nu_j} - \eta ( 1 -\sum_{dx_i=1} \mu_i ) - \delta ( 1 -\sum_{dy_j=1} \nu_j ) - \lambda ( \mu^T H_1 \nu, \ldots , \mu^T H_p \nu )^T $$ where the variables are defined as follows:

\(x\) is a vector of uncensored data for the first sample

\(y\) is a vector of uncensored data for the second sample

\(wx\) is a vector of estimated weights for the first sample

\(wy\) is a vector of estimated weights for the second sample

\(\mu\) is a vector of estimated probability jumps for the first sample

\(\nu\) is a vector of estimated probability jumps for the second sample

\(H_k = [ g_k(x_i, y_j) - mean_k ], k=1, \ldots , p\), where \(g_k(x,y)\) is a user-chosen function

\(H = [H_1, ... , H_p]\) (used as argument in el.cen.EMm function, which calls el2.test.wtm)

\(mean\) is a vector of length \(p\) of hypothesized means, such that \(mean_k\) is the hypothesized value of \(E(g_k(x,y))\)

\(E\) indicates ``expected value''

Usage

el2.test.wtm(xd1,yd1,wxd1new, wyd1new, muvec, nuvec, Hu, Hmu, Hnu, p, mean, maxit=15)

Value

el2.test.wtm returns a list of values as follows:

constmat

a matrix of row vectors, where the \(k\)th row vector holds successive values of \(\mu^T H_k \nu , k=1, \ldots, p\), generated by the Newton-Raphson algorith m

lam

the vector of Lagrangian mulipliers used in the calculations

muvec1

the vector of probability jumps for the first sample that maximize the weighted empirical likelihood

nuvec1

the vector of probability jumps for the second sample that maximize the weighted empirical likelihood

mean

the vector of hypothesized means

Arguments

xd1

a vector of uncensored data for the first sample

yd1

a vector of uncensored data for the second sample

wxd1new

a vector of estimated weights for xd1

wyd1new

a vector of estimated weights for yd1

muvec

a vector of estimated probability jumps for xd1

nuvec

a vector of estimated probability jumps for yd1

Hu

\(H_u = [ H_1 - [mean_1], \ldots , H_p - [mean_p] ], dx_i=1, dy_j=1\)

Hmu

a matrix, whose calculation is shown in the example below

Hnu

a matrix, whose calculation is shown in the example below

p

the number of hypotheses

mean

a vector of hypothesized values of \(E(g_k(u,v)), k=1, \ldots,p\)

maxit

a positive integer used to control the maximum number of iterations in the Newton-Raphson algorithm; default is 10

Author

William H. Barton <[e-mail]>

Details

This function is called by el2.cen.EMm. It is listed here because the user may find it useful elsewhere.

The value of \(mean_k\) should be chosen between the maximum and minimum values of \(g_k(xd1_i,yd1_j)\); otherwise there may be no distributions for \(xd1\) and \(yd1\) that will satisfy the the mean-type hypothesis. If \(mean_k\) is inside this interval, but the convergence is still not satisfactory, then the value of \(mean_k\) should be moved closer to the NPMLE for \(E(g(xd1,yd1))\). (The NPMLE itself should always be a feasible value for \(mean_k\).) The calculations for this function are derived in Owen (2001).

References

Owen, A.B. (2001). Empirical Likelihood. Chapman and Hall/CRC, Boca Raton, pp.223-227.

Examples

#Ho1: P(X>Y) = 0.5
#Ho2: P(X>1060) = 0.5
#g1(x) = I[x > y]
#g2(y) = I[x > 1060]

mean<-c(0.5,0.5)
p<-2

xd1<-c(10,85,209,273,279,324,391,566,852,881,895,954,1101,1393,1669,1823,1941)
nx1=length(xd1)
yd1<-c(21,38,39,51,77,185,524,610,612,677,798,899,946,1010,1074,1147,1154,1329,1484,1602,1952)
ny1=length(yd1)

wxd1new<-c(2.267983, 1.123600, 1.121683, 1.121683, 1.121683, 1.121683, 1.121683,
 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.261740, 2.912753,
 2.912753, 2.912753)
muvec<-c(0.08835785, 0.04075290, 0.04012084, 0.04012084, 0.04012084, 0.04012084,
 0.04012084, 0.03538020, 0.03389263, 0.03389263, 0.03389263, 0.03322693, 0.04901516,
 0.05902008, 0.13065491, 0.13065491, 0.13065491)

wyd1new<-c(1.431653, 1.431653, 1.431653, 1.431653, 1.431653, 1.438453, 1.079955, 1.080832,
 1.080832, 1.080832, 1.080832, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000,
 1.222883, 1.227865, 1.739636, 5.809616)
nuvec<-c(0.04249966, 0.04249966, 0.04249966, 0.04249966, 0.04249966, 0.04316922, 0.03425722,
 0.03463312, 0.03463312, 0.03463312, 0.03463312, 0.03300598, 0.03300598, 0.03333333,
 0.03333333, 0.03382827, 0.03382827, 0.04136800, 0.04229270, 0.05992020, 0.22762676)


H1u<-matrix(NA,nrow=nx1,ncol=ny1)
H2u<-matrix(NA,nrow=nx1,ncol=ny1)
for (i in 1:nx1) {
   for (j in 1:ny1) {
        H1u[i,j]<-(xd1[i]>yd1[j])
        H2u[i,j]<-(xd1[i]>1060) } }
Hu=matrix(c(H1u,H2u),nrow=nx1,ncol=p*ny1)
for (k in 1:p) {
     M1 <- matrix(mean[k], nrow=nx1, ncol=ny1)
     Hu[,((k-1)*ny1+1):(k*ny1)] <- Hu[,((k-1)*ny1+1):(k*ny1)] - M1}
Hmu <- matrix(NA,nrow=p, ncol=ny1*nx1)
Hnu <- matrix(NA,nrow=p, ncol=ny1*nx1) 
for (i in 1:p) {
   for (k in 1:nx1) {
        Hmu[i, ((k-1)*ny1+1):(k*ny1)] <-  Hu[k,((i-1)*ny1+1):(i*ny1)] } }
for (i in 1:p) {
   for (k in 1:ny1) {
        Hnu[i,((k-1)*nx1+1):(k*nx1)] <- Hu[(1:nx1),(i-1)*ny1+k]} }

el2.test.wtm(xd1,yd1,wxd1new, wyd1new, muvec, nuvec, Hu, Hmu,
  Hnu, p, mean, maxit=10)

#muvec1
# [1] 0.08835789 0.04075290 0.04012083 0.04012083 0.04012083 0.04012083 0.04012083
# [8] 0.03538021 0.03389264 0.03389264 0.03389264 0.03322693 0.04901513 0.05902002
# [15] 0.13065495 0.13065495 0.13065495

#nuvec1
# [1] 0.04249967 0.04249967 0.04249967 0.04249967 0.04249967 0.04316920 0.03425722
# [8] 0.03463310 0.03463310 0.03463310 0.03463310 0.03300597 0.03300597 0.03333333
# [15] 0.03333333 0.03382827 0.03382827 0.04136801 0.04229269 0.05992018 0.22762677

#  $lam
#        [,1]      [,2]
#  [1,] 8.9549 -10.29119