Re: [R] Improving effeciency - better table()?

From: Simon Cullen <cullens_at_tcd.ie>
Date: Wed 07 Jul 2004 - 02:21:44 EST

It seems that attachments are stripped off so I'll include the code and the output of Rprof below:

#apologies for disgusting wrapping

f.powertest <-
function(path=paste(substr(getwd(),1,3),"R/",sep=""),iter=5000,sample=25,trun.prop=0,m,s){

        rej.2.5 <- 0                         

	for (i in 1:iter){
		z<-c()
		z.trun<-c()
		zeros<- 0
		total<-0
		while(length(z.trun)< sample){
			
			trun<-qnorm(trun.prop,m,s)
			z<-rnorm(sample,m,s)
			ind <- z>trun
			zeros <- zeros + sum(!ind)
			total<- total + sample
			z.trun <- c(z.trun,z[ind])
			
		}

		z.trun<-z.trun[1:sample]
		tau<-ifelse(qnorm(zeros/total)==-Inf,-3,qnorm(zeros/total))
		
		#estimate the true mean
		m.tr <- mean(z.trun)
		s.tr <- sd(z.trun)
		s.untr <- s.tr/sqrt(1-f.lambda(tau)*(f.lambda(tau)-tau))
		m.untr <- m.tr - s.untr*f.lambda(tau)

		#cacluate the breakpoints
		
		#number of breakpoints first:
		br <- ceiling(log2(length(z.trun)) + 1)
		z.trun<- sort(z.trun)
	
		br.len <- (z.trun[sample - 6] - z.trun[1])/(br-1)
		
		breaks <- seq(from=z.trun[1],by=br.len,length=br)
		top <- 1000#top limit of the range - I'm using N(0,1) so for the moment  
this is ok
		breaks <- c(breaks, top)
		tab <- hist(z.trun,br=breaks, plot=F,right=F)$counts
		expect<-f.expect(breaks,m.untr,s.untr,tau)*sample

		stat <- sum((tab-expect)^2/expect)

		if(qchisq(0.975,br-3) < stat){
			rej.2.5 <-rej.2.5 + 1
		}

	}

	return(100*rej.2.5/iter)

}

f.expect <- function(breaks,m,sigma,alpha){
	res<-sapply(breaks,f.tr.pnorm,m=m,sigma=sigma,alpha=alpha)
	return(res[-1] - res[-length(res)])

}

f.tr.dnorm <- function(x,m=0,sigma=1,alpha){

        dnorm((x-m)/sigma)/(sigma*(1-pnorm(alpha)))
}

f.tr.pnorm <- function(q,m=0,sigma=1,alpha){

        integrate(f.tr.dnorm,alpha,q,m=m,sigma=sigma,alpha=alpha)$value#rel.tol=.Machine$double.eps^0.5
}

f.test <- function(){

	Rprof()
	f.powertest(iter=1000,m=0,s=1,sample=200,trun.prop=0.25)
	Rprof(NULL)
	print(summaryRprof())

}

f.lambda <- function(x) dnorm(x)/(1 - pnorm(x))


> f.test()
$by.self

                    self.time self.pct total.time total.pct
integrate               0.92     15.2       3.28      51.7
as.double               0.56      9.3       0.68      10.7
dnorm                   0.44      7.3       0.52       8.2
rnorm                   0.30      5.0       0.30       4.7
pnorm                   0.26      4.3       0.26       4.1
>                       0.18      3.0       0.18       2.8
hist.default            0.18      3.0       1.24      19.6
f                       0.16      2.6       0.90      14.2
names                   0.16      2.6       0.18       2.8
f.powertest             0.14      2.3       6.34     100.0
switch                  0.14      2.3       0.14       2.2
==                      0.12      2.0       0.12       1.9
as.double.default       0.12      2.0       0.12       1.9
as.integer              0.12      2.0       0.16       2.5
diff.default            0.12      2.0       0.32       5.0
-                       0.10      1.7       0.10       1.6
paste                   0.10      1.7       0.20       3.2
seq                     0.10      1.7       0.20       3.2
...

$by.total

                    total.time total.pct self.time self.pct
f.powertest              6.34     100.0      0.14      2.3
f.test                   6.34     100.0      0.00      0.0
f.expect                 3.64      57.4      0.02      0.3
sapply                   3.62      57.1      0.04      0.7
lapply                   3.52      55.5      0.04      0.7
FUN                      3.38      53.3      0.02      0.3
integrate                3.28      51.7      0.92     15.2
hist                     1.30      20.5      0.06      1.0
hist.default             1.24      19.6      0.18      3.0
<Anonymous>              0.94      14.8      0.04      0.7
f                        0.90      14.2      0.16      2.6
as.double                0.68      10.7      0.56      9.3
dnorm                    0.52       8.2      0.44      7.3
sort                     0.40       6.3      0.04      0.7
diff                     0.36       5.7      0.00      0.0
storage.mode<-           0.34       5.4      0.02      0.3
diff.default             0.32       5.0      0.12      2.0
rnorm                    0.30       4.7      0.30      5.0
pnorm                    0.26       4.1      0.26      4.3
ifelse                   0.22       3.5      0.04      0.7
eval                     0.20       3.2      0.06      1.0
paste                    0.20       3.2      0.10      1.7
seq                      0.20       3.2      0.10      1.7
...
$sampling.time
[1] 6.34
-- 
SC

Simon Cullen
Room 3030
Dept. Of Economics
Trinity College Dublin

Ph. (608)3477
Email cullens@tcd.ie

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Received on Wed Jul 07 02:22:22 2004

This archive was generated by hypermail 2.1.8 : Wed 03 Nov 2004 - 22:54:44 EST