[R] Function to approximate complex integral

[Douglas Bates]

On 4/19/06, Doran, Harold <HDoran at air.org> wrote:

> I am writing a small function to approximate an integral that cannot be
> evaluated in closed form. I am partially successful at this point and am
> experiencing one small, albeit important problem. Here is part of my
> function below.

> This is a psychometric problem for dichotomously scored test items where
> x is a vector of 1s or 0s denoting whether the respondent answered the
> item correctly (1) or otherwise (0), b is a vector of item difficulties,
> and theta is an ability estimate for the individual.
>
> rasch <- function(b,theta){
>    1 / ( 1 + exp(b - theta))
>    }
>
> The function rasch gives the probability of a correct response to item i
> conditional on theta, the individuals ability estimate
>
> myfun <- function(x, b, theta){
>    sum(rasch(b, theta)^x * ( 1 - rasch(b,theta) )^(1-x) * dnorm(theta))
>    }
>
> This is the likelihood function assuming the data are Bernoulli
> distributed multiplied by a population distribution.
>
> Now, when x,b, and theta are of equal length the function works fine as
> below
> x <- c(1,1,0)
> b <-   c(-2,-1,0)
> t <-   c(-2,-1.5,-1)
> > myfun(x,b,t)
> [1] 0.2527884
>
> However, I want theta to be a vector of discrete values that will be
> larger than both x and b, something like
>
> t <- seq(-5, 0, by = .01)
>
> However, this gives me the error
> > myfun(x,b,t)
>
>  Warning messages:
> 1: longer object length
>         is not a multiple of shorter object length in: b - theta
>
> So, for the problem above, I want item 1 to be evaluated at theta 1
> through theta q and then item 2 is evaluated at theta 1 and through
> theta q and so forth for each item.
>
> Can anyone recommend a way for me to modify my function above to
> accomplish this?

I believe (partially based on our off-list discussions of this model)
that you want to use the function "outer", as in

rasch <- function(b, theta) 1/(1 - exp(outer(theta, b, "-")))

If length(theta) == m and length(b) == n then this returns the m by n
matrix of probabilities of subject i correctly answering question j.