Statalist


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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


From   "Martin Weiss" <martin.weiss1@gmx.de>
To   <statalist@hsphsun2.harvard.edu>
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)
#   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


bdwamena@umich.edu

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" <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/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index