<>
I am forwarding this message on behalf of Ben.
HTH
Martin
_______________________
----- Original Message -----
From: "Ben Dwamena" <[email protected]>
To: "Martin Weiss" <[email protected]>
Sent: Saturday, May 30, 2009 2:11 PM
Subject: st: Re: Translating an R expression (diff <- pchisq(d,
p) -(0.5:n)/n) into Stata code
Hi Martin,
This is the full code I was translating:
function(x,m0,c0,alpha,pcrit){
# Adaptive reweighted estimator for multivariate location and scatter
# with hard-rejection weights and delta = chi2inv(1-d,p)
#
# Input arguments
# x: Dataset (n x p)
# m0: Initial location estimator (1 x p)
# c0: Initial scatter estimator (p x p)
# alpha: Maximum thresholding proportion
# (optional scalar, default: alpha = 0.025)
# pcrit: critical value for outlier probability
# (optional scalar, default values from simulations)
#
# Output arguments:
# m: Adaptive location estimator (p x 1)
# c: Adaptive scatter estimator (p x p)
# cn: Adaptive threshold (scalar)
# w: Weight vector (n x 1)
#
n <- nrow(x)
p <- ncol(x)
# Critical value for outlier probability based on simulations for
alpha=0.025
if (missing(pcrit)){
if (p<=10) pcrit <- (0.24-0.003*p)/sqrt(n)
if (p>10) pcrit <- (0.252-0.0018*p)/sqrt(n)
}
if (missing(alpha)) delta<-qchisq(0.975,p) else delta<-qchisq(1-alpha,p)
d2<-mahalanobis(x,m0,c0)
d2ord <- sort(d2)
dif <- pchisq(d2ord,p) - (0.5:n)/n
i <- (d2ord>=delta) & (dif>0)
if (sum(i)==0) alfan<-0 else alfan<-max(dif[i])
if (alfan<pcrit) alfan<-0
#if (alfan>0) cn<-max(d2ord[n-floor(n*alfan)],delta) else cn<-Inf
if (alfan>0) cn<-max(d2ord[n-ceiling(n*alfan)],delta) else cn<-Inf
w <- d2<cn
if(sum(w)==0) {
m <- m0
c <- c0
}
else {
m <- apply(x[w,],2,mean)
c1 <- as.matrix(x-rep(1,n)%*%t(m))
c <- (t(c1)%*%diag(w)%*%c1)/sum(w)
}
list(m=m,c=c,cn=cn,w=w)
}
Thanks,
Ben
Ben Adarkwa Dwamena, MD
Assistant Professor of Radiology
Division of Nuclear Medicine
Department of Radiology
University of Michigan Health Systems
B1G505 UH SPC 5028
1500 E. Medical Center Drive
Ann Arbor, MI 48109-5028
[email protected]
http://sitemaker.umich.edu/metadiagnosis
Staff Physician
D748 Nuclear Medicine Service (115),
VA Ann Arbor Health Care System
2215 Fuller Road
Ann Arbor, MI 48105
734-845-3775 Phone
734-845-3252 Fax
"Martin Weiss" <[email protected]> 5/30/2009 4:31 AM >>>
<>
Here is an attempt to translate your code. What is "d", btw?
***
clear*
local d 4
//# of rows
set obs 100
//# of columns set to 3
gen x=.
gen y=.
gen z=.
forv i=0.05(0.1)`=c(N)/10'{
di chi2(`d', c(k))-`i'
}
***
HTH
Martin
_______________________
----- Original Message -----
From: "Ben Dwamena" <[email protected]>
To: <[email protected]>
Sent: Saturday, May 30, 2009 8:28 AM
Subject: st: Translating an R expression (diff <- pchisq(d, p) - (0.5:n)/n)
into Stata code
Dear Listers with R experience,
What is an equivalent Stata code for the following R expression:
diff <- pchisq(d, p) - (0.5:n)/n
where x is dataset (n x p)
n is number of rows of x
p is number of columns of x
Thanks,
Ben
**********************************************************
Electronic Mail is not secure, may not be read every day, and should not
be used for urgent or sensitive issues
*
* For searches and help try:
* http://www.stata.com/help.cgi?search
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
*
* For searches and help try:
* http://www.stata.com/help.cgi?search
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
**********************************************************
Electronic Mail is not secure, may not be read every day, and should not be
used for urgent or sensitive issues
*
* For searches and help try:
* http://www.stata.com/help.cgi?search
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
