# st: Fw: Re: Translating an R expression (diff <- pchisq(d, p) -(0.5:n)/n) into Stata code

 From "Martin Weiss" To Subject st: Fw: Re: Translating an R expression (diff <- pchisq(d, p) -(0.5:n)/n) into Stata code Date Sat, 30 May 2009 15:13:31 +0200

<>

I am forwarding this message on behalf of Ben.

HTH
Martin
_______________________
----- Original Message ----- From: "Ben Dwamena" <bdwamena@med.umich.edu>
To: "Martin Weiss" <martin.weiss1@gmx.de>
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)
#   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

Division of Nuclear Medicine
University of Michigan Health Systems
B1G505 UH SPC 5028
1500 E. Medical Center Drive
Ann Arbor, MI  48109-5028

bdwamena@umich.edu

Staff Physician
D748 Nuclear Medicine Service (115),
VA Ann Arbor Health Care System
Ann Arbor, MI 48105
734-845-3775 Phone
734-845-3252 Fax

"Martin Weiss" <martin.weiss1@gmx.de> 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" <bdwamena@med.umich.edu>
To: <statalist@hsphsun2.harvard.edu>
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/