el2.test.wtm {emplik2} | R Documentation |
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{μ_i} + sum_{dy_j=1} wy_j log{nu_j} - eta ( 1 -sum_{dx_i=1} μ_i ) - delta ( 1 -sum_{dy_j=1} nu_j ) - λ ( μ^T H_1 nu, ... , μ^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
μ 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, ... , 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''
el2.test.wtm(xd1,yd1,wxd1new, wyd1new, muvec, nuvec, Hu, Hmu, Hnu, p, mean, maxit=10)
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], ... , 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, ...,p |
maxit |
a positive integer used to control the maximum number of iterations in the Newton-Raphson algorithm; default is 10 |
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).
el2.test.wtm
returns a list of values as follows:
constmat |
a matrix of row vectors, where the kth row vector holds successive values of μ^T H_k nu , k=1, ..., 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 |
William H. Barton <bbarton@lexmark.com>
Owen, A.B. (2001). Empirical Likelihood
. Chapman and Hall/CRC, Boca Raton, pp.223-227.
#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