I am trying to calculate the marginal maximum likelihood estimate of a parameter ρ given a vector of data x for person i across n items. This is in the context of Item Response Theory. The first (usual IRT) parameter, θ, is the ability (theoretically unbounded but usually between -3 and 3). The second parameter, ρ, is the probability that the respondent answers in a non-optimal way. I want to integrate θ out of my joint likelihood equation:
where the probability function in the integrand is
The second term in the integrand, f(θ), will just be a numerical approximation to the probability density function of θ. The three dots in the first set of curly brackets above can take three different forms depending on certain conditions: 0 (in the "normal" condition, so person i will definitely not respond in a non-optimal way), the formula in the second set of curly brackets, or one minus the formula in the second set of curly brackets.
My problem is just in wrapping my head around how I can accomplish these calculations. I am using R, but I suppose any good explanation of how I can do this computationally would be sufficient for me to understand enough to try to implement it. It is my understanding that in order to be able to carry out these calculations, I need to express the first term in the integrand (i.e., the joint probability function) as a function of 0. Then, I multiply this function of 0 by many numerical approximations to the second term (pdf of 0). Finally, I integrate each of these products over 0. Does this mean that the integrand is entirely a function of 0, and once I integrate 0 out I essentially go from a formula to a number? At which point I would sum of the results over all n items, and then implement an optimization routine to search for the most likely value of ρ within the bounds allowed?