From: Doran, Harold <HDoran_at_air.org>

Date: Sat 06 May 2006 - 01:54:55 EST

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Sun May 07 06:26:29 2006

Date: Sat 06 May 2006 - 01:54:55 EST

Spencer

Thanks for your thoughts on this. I did a bit of work and did end up with a method (more a trick), but it did work. I am certain there are better ways to do this, but here is how I resolved the issue.

The integral I need to evaluate is

\begin{equation}

\frac{\int_c^{\infty} p(x|\theta)f(\theta)d\theta}
{\int_{-\infty}^{\infty} p(x|\theta)f(\theta)d\theta}
\end{equation}

Where p(x|\theta) is the likelihood for the Rasch item response theory
model and f(\theta) is the population

distribution. Originally, I was using numerical quadrature by summing
from -10 to 10 in increments of .01 for the denominator and from c to 10
for the numerator. This worked, but was rather expensive.

Now, the gauss.quad.prob function provides weights that I can use for
the denominator via:

> gauss.quad.prob(49, dist='normal', mu=mu, sigma=sigma)

Where mu and sigma are values that are estimated from another procedure. But, for sake of argument assume N(0,1). My challenge was to get the appropriate nodes and weights for the numerator since I am going from -Inf to c. So, what I did to get this to work was get the quadrature nodes and weights using

> gauss.quad(49,kind="laguerre")

Then, I do the following to get them to work with a normal distribution:

F(y) = p(y+c) * N(y+c|\mu, \sigma^2)

Where y is the laguerre node and c is the lower bound of integration. Subsequently, I plug this into the following

\sum_{i=1}^Q F(y_i)e^y_i \cdot w_i

Where i indexes node and w is the weight at node y_i.

As it turns out, I can evaluate this using only 49 quadrature points (and not the 2001 I was using with equal intervals) and get my result. Both provide the same answer, but CPU time is reduced to basically zero for the quadrature whereas it was about 4 CPU seconds for the equal interval approach. Keep in mind that I am looping through a student achievement data file where this is being applied to about 70,000 students, so time matters here.

Here is how I programmed this. I realize there are much better ways to accomplish this, and any suggestions would be very much appreciated.

library(statmod)

rasch <- function(b,theta){

1 / (1 + exp(outer(b,theta,'-')))

}

like.mat <- function(x,b,theta){

rasch(b, theta)^x * (1 - rasch(b,theta))^(1-x) }

class.numer <- function(x,b, prof_cut, mu=0, sigma=1, aboveQ){

gauss_numer <- gauss.quad(49,kind="laguerre") if(aboveQ==FALSE){

mat <- rbind(like.mat(x,b, (prof_cut-gauss_numer$nodes)), dnorm(prof_cut-gauss_numer$nodes, mean=mu, sd=sigma))

} else { mat <- rbind(like.mat(x,b, (gauss_numer$nodes+prof_cut)), dnorm(gauss_numer$nodes+prof_cut, mean=mu, sd=sigma))

}

f_y <- rbind(apply(mat, 2, prod), exp(gauss_numer$nodes),
gauss_numer$weights)

sum(apply(f_y,2,prod))

}

class.denom <- function(x,b, mu=0, sigma=1){

gauss_denom <- gauss.quad.prob(49, dist='normal', mu=mu, sigma=sigma)
mat <- rbind(like.mat(x,b,gauss_denom$nodes),gauss_denom$weights)
sum(apply(mat, 2, prod))

}

class.acc <-function(x,b,prof_cut, mu=0, sigma=1, aboveQ=TRUE){

result <- class.numer(x,b,prof_cut, mu,sigma, aboveQ)/class.denom(x,b, mu, sigma)

return(result)

}

# test the function

set.seed(1)

response <- c(rep(1,10), rep(0,10))

items <- runif(20, -3,3)

class.acc(response, items, .05, above=F)
[1] 0.3279967

**> GAUSSIAN QUADRATURE?
**

>

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Sun May 07 06:26:29 2006

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.1.8, at Sun 07 May 2006 - 12:10:01 EST.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*