similar to: Basic vector operations was: Function to approximate complex integral

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

This is a very basic question, but I am a bit confused with optim. I want to get the MLEs using optim which could replace the newton-raphson code I have below which also gives the MLEs. The function takes as input a vector x denoting whether a respondent answered an item correctly (x=1) or not (x=0). It also takes as input a vector b_vector, and these are parameters of test items (Rasch estimates

I have constructed the function mml2 (below) based on the likelihood function described in the minimal latex I have pasted below for anyone who wants to look at it. This function finds parameter estimates for a basic Rasch (IRT) model. Using the function without the gradient, using either nlminb or optim returns the correct parameter estimates and, in the case of optim, the correct standard

I have a function that I am currently using very inefficiently. The following are needed to illustrate the problem: set.seed(12345) dat <- matrix(sample(c(0,1), 110, replace = TRUE), nrow = 11, ncol=10) mis <- sample(1:110, 5) dat[mis] <- NA theta <- rnorm(11) b_vector <- runif(10, -4,4) empty <- which(is.na(t(dat))) So, I have a matrix (dat) with some values within the matrix

List: Thank you for the replies to my post yesterday. Gabor and Phil also gave useful replies on how to improve the function by relying on mapply rather than the explicit for loop. In general, I try and use the family of apply functions rather than the looping constructs such as for, while etc as a matter of practice. However, it seems the mapply function in this case is slower (in terms of CPU

Dear list Suppose I have the following vector: x <- c(3,4,2,5,6) Obviously, this sums to 20. Now, I want to have a second vector, call it x2, that sums to x where 5 <= x <= 20, but there are constraints. 1) The new vector must be same length as x 2) No element of the new vector can be 0 3) Element x2[i] of the new vector cannot be larger than element x[i] of the original vector 4)

I need to wrap a loop inside a function and am having a small bit of difficulty getting the results I need. Below is a replicable example. # define functions pcm <- function(theta,d,score){ exp(rowSums(outer(theta,d[1:score],'-')))/ apply(exp(apply(outer(theta,d, '-'), 1, cumsum)), 2, sum) } foo <- function(theta,items, score){ like.mat <-

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

I have the following item response theory (IRT) likelihood that I want to maximize w.r.t. to theta (student ability). L(\theta) = \prod(p(x)) Where p(x) is the 3-parameter logistic model when items are scored dichotomously (x_{ij} = 0 or 1) and p(x) is Muraki's generalized partial credit model when items are scored polytomously (x_{ij} = 0 \ldots J). Now, I wrote the following two functions

Suppose I have a function, such as the toy example below: myFun <- function(x, max.iter = 5) { for(i in 1:10){ result <- x + i iter <- i if(iter == max.iter) stop('Max reached') } result } I can of course do this: myFun(10, max.iter = 11) However, if I reach the maximum

Dear List: I have created some code to simulate data from a complex sample where 5000 students are nested in 50 schools. My code returns a dataframe with a variable representing student achievement at a single time point. My actual code for creating this is below. What I would like to do is generate a second column of data that is correlated with the first at .8 and has the same means within

Dear List: We are generating data such that students are clustered in schools for some item response data for a simulation study. One component of our simulation is to generate measurement error from a logistic distribution with a mean of 0 and standard deviation of 1.7 to match the logistic curve of the Rasch model. We are generating an error term for each of the 40 hypothetical test items a

Hello R, I am looking for non-parametric simulation in IRT. Is there any IRT package that does non-parametric simulation? helen L [[alternative HTML version deleted]]

Suppose I have the following function myFun <- function(formula, data){ f <- formula(formula) dat <- model.frame(f, data) dat } Applying it with this sample data yields a new dataframe: qqq <- data.frame(grade = c(3, NA, 3,4,5,5,4,3), score = rnorm(8), idVar = c(1:8)) dat <- myFun(score ~ grade, qqq) However, what I would like is for the resulting dataframe (dat) to include

I am using both nlminb and optim to get MLEs from a likelihood function I have developed. AFAIK, the model I has not been previously used in this way and so I am struggling a bit to unit test my code since I don't have another data set to compare this kind of estimation to. The likelihood I have is (in tex below) \begin{equation} \label{eqn:marginal} L(\beta) = \prod_{s=1}^N \int

I have just submitted MiscPsycho to CRAN. MiscPsycho contains functions for miscellaneous psychometrics that may be useful for applied psychometricians. MML estimation already exists in the ltm package. Hence, a jml option is provided for users who prefer this method. The jml function gives back rasch difficulties and the same Infit and Outfit statistics as Winsteps. Also, jml is known to return

I have just submitted MiscPsycho to CRAN. MiscPsycho contains functions for miscellaneous psychometrics that may be useful for applied psychometricians. MML estimation already exists in the ltm package. Hence, a jml option is provided for users who prefer this method. The jml function gives back rasch difficulties and the same Infit and Outfit statistics as Winsteps. Also, jml is known to return

Does anyone see why my code does not integrate to 1? library(statmod) mu <- 0 s <- 1 Q <- 5 qq <- gauss.quad(Q, kind='hermite') sum((1/(s*sqrt(2*pi))) * exp(-((qq$nodes-mu)^2/(2*s^2))) * qq$weights) ### This does what's it is supposed to myNorm <- function(theta) (1/(s*sqrt(2*pi))) * exp(-((theta-mu)^2/(2*s^2))) integrate(myNorm, -Inf, Inf)

Dear list: I'm hoping to tap in to the statistical expertise in the group, especially those familiar with simulation techniques. I'm finalizing a study where I obtain standard errors from two sources. The first source is a monte carlo simulation and the other source is an analytical model I have developed that appears to recover the standard errors from the simulation. All analysis are

Dear friends, I'm writing a function with three arguments myfun <- function(theta, X, values) { .... .... } in this function, I'm trying to write consistency checks. In order to compute the statistic of interest I only need theta and values. The idea of having X in there is that, if values is not provided by the user, then values is computed from X. my problem is I'm trying to