求助解决自编参数IRT联合极大似然估计错误

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

我编了一个单参数IRT的联合极大释然估计函数JMLE.1P（），但是报错：Error in JMLE.1P(data1$U_Matrix) : object 'theta.temp' not found  。其中data1$U_Matrix是作答矩阵。请问我是哪里出错了呢？谢谢！求大家指点！其中CMLE.theta.1P（）和CMLE.beta.1P()这两个函数是我已经编好的能力参数和难度参数的条件极大似然估计函数。JMLE.1P()会调用上面两个函数。

1. #1PLM/Rasch moedel JMLE
   
2. #作答矩阵需事先排除满分和零分的被试数据，以及全对和全错的项目数据
   
3. 
   
4. if (!exists("CMLE.theta.1P")) source("CMLE.theta.1P.R")
   
5. if (!exists("CMLE.beta.1P"))  source("CMLE.beta.1P.R")
   
6. 
   
7. JMLE.1P<- function(Resp.Matrix,n2=20,tol2=0.005){
   
8.   N<- nrow(Resp.Matrix)
   
9.   m<- ncol(Resp.Matrix)
   
10.   T<- rowSums(Resp.Matrix)
    
11.   mean_T<- mean(T)
    
12.   sd_T<- sd(T)
    
13.   p<- numeric()
    
14.   for (j in 1:m) p[j]<- sum(Resp.Matrix[,j])/N
    
15.   q<- 1-p
    
16.   
    
17.   #取在题目上得分的被试的卷面总分平均分，一共有m道题，mean_Tp就有m个元素
    
18.   mean_Tp<- numeric()
    
19.   for (j in 1:m) mean_Tp[j]<- sum(Resp.Matrix[,j]*T)/sum(Resp.Matrix[,j])
    
20.   
    
21.   #在题目上失分的被试的卷面总分平均分
    
22.   mean_Tq=numeric()
    
23.   for (j in 1:m) mean_Tq[j]<- sum((1-Resp.Matrix[,j])*T)/sum(1-Resp.Matrix[,j])
    
24.   
    
25.   #设定beta初始值
    
26.   #有题目通过率极端计算对应的概率密度，用于将后面的点二列相关转换为二列相关，进而设定项目参数估计初值
    
27.   y<- dnorm(qnorm(p,lower.tail=F))
    
28.   r<- (mean_Tp-mean_Tq)*(p*q)/(sd_T*y) #二列相关
    
29.   beta.start<- qnorm(p,lower.tail=F)/r
    
30.   beta.temp<- beta.start
    
31.   theta.start<- log(T/(m-T))
    
32.   P.start<- matrix(,N,m)
    
33.   for (i in 1:N) {
    
34.     for(j in 1:m) {
    
35.       P.start[i,j]<- 1/(1 + exp(-1.7 * (theta.start[i] - beta.start[j])))
    
36.     }
    
37.   }
    
38.   logLik.start<-  sum(Resp.Matrix*log(P.start)+(1-Resp.Matrix)*log(1-P.start))
    
39. 
    
40.   
    
41.   logLikDiff<- 999999999
    
42.   iterN<- 1
    
43.   while (logLikDiff < tol2 & iterN < n2){
    
44.     cmle.theta.temp<- apply(Resp.Matrix,1,CMLE.theta.1P,beta=beta.temp)
    
45.     theta.temp<- numeric()
    
46.     for (i in 1:N) theta.temp[i]<- cmle.theta.temp[[i]]> JMLE.1P(data1$U_Matrix)
    
47. Error in JMLE.1P(data1$U_Matrix) : object 'theta.temp' not found
    
48. point estimator'
    
49.    
    
50.     cmle.beta.temp<- apply(Resp.Matrix,2,CMLE.beta.1P,theta=theta.temp,k=20)
    
51.     beta.temp<- numeric()
    
52.     for(j in 1:m) beta.temp[j]<- cmle.beta.temp[[j]]> JMLE.1P(data1$U_Matrix)
    
53. Error in JMLE.1P(data1$U_Matrix) : object 'theta.temp' not found
    
54. point estimator'
    
55.    
    
56.     P.temp<- matrix(,N,m)
    
57.     for (i in 1:N) {
    
58.       for(j in 1:m) {
    
59.         P.temp[i,j]<- 1/(1 + exp(-1.7 * (theta.temp[i] - beta.temp[j])))
    
60.       }
    
61.     }
    
62.     logLik.temp<- sum(Resp.Matrix*log(P.temp)+(1-Resp.Matrix)*log(1-P.temp))
    
63.     logLikDiff<- abs(logLik.start-logLik.temp)
    
64.     iterN<- iterN+1
    
65.   }
    
66.   Results<- list(theta.estimator=theta.temp, beta.estimator=beta.temp)
    
67.   return(Results)
    
68. }

复制代码

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：极大似然估计 似然估计 极大似然 IRT matrix source Error 联合 能力 项目

超级大菜鸟 发表于 2014-11-27 17:14
我编了一个单参数IRT的联合极大释然估计函数JMLE.1P（），但是报错：Error in JMLE.1P(data1$U_Matrix) : o ...

自己顶，各位大神帮帮忙啊！！！！

Resp.Matrix 和 data1$U_Matrix 是同样的吗？

CMLE.theta.1P（）和CMLE.beta.1P(） 这两个函数不在，所以很难运行你的Code 帮你查问题出在哪里

calsunny 发表于 2014-12-2 14:29
Resp.Matrix 和 data1$U_Matrix 是同样的吗？

CMLE.theta.1P（）和CMLE.beta.1P(） 这两个函数不在，所 ...

1. CMLE.beta.1P<- function(resp.vect,D=1.7,theta,beta.start=0,k=20,n1=50,tol1=0.0005,...){
   
2.   beta<- beta.start
   
3.   d1.logLikDiff<- 99999999999
   
4.   it.log<- matrix(ncol=3,nrow=0)
   
5.   colnames(it.log)<- c("beta_estimator","d1.logLik","d2.logLik")
   
6.   
   
7.   #将某题的作答向量与theta合并成一个数据框
   
8.   A<- data.frame(resp.vect,theta)
   
9.   A<- A[order(A$theta),] #以A中的theta排序
   
10.   A$theta_space<- cut(A$theta,breaks=k)
    
11.   
    
12.   #将theta分成k组后，取每组的中位数代表该组的能力水平
    
13.   med_theta_space<- aggregate(A$theta,list(A$theta_space),median)
    
14.   names(med_theta_space)<- c("theta_space","theta_level")
    
15.   
    
16.   #B是每个能力水平答对r和答错(f-r)组成的矩阵
    
17.   B<- table(A$resp.vect,A$theta_space)
    
18.   f<- colSums(B)
    
19.   r<- B[2,]
    
20.   
    
21.   iterN<- 1 #iteration number
    
22.   
    
23.   while(d1.logLikDiff > tol1 & iterN< n1) {
    
24.     probs.1<- 1/(1+exp(-D*(med_theta_space$theta_level - beta)))
    
25.     probs.0<- 1-probs.1
    
26.    
    
27.     #First derivarive
    
28.     d1.logLik<- sum(D*(f*probs.1-r))
    
29.    
    
30.     #Second derivative
    
31.     d2.logLik<- -sum(D*f*probs.1*probs.0)
    
32.     it.log<- rbind(it.log, c(beta, d1.logLik, d2.logLik))
    
33.    
    
34.     #Newton Raphson
    
35.     beta<- beta - (d1.logLik/d2.logLik)
    
36.    
    
37.     #Convergence  critieria   
    
38.     if(nrow(it.log) > 1) {
    
39.       d1.logLikDiff<- abs(d1.logLik-it.log[(nrow(it.log) - 1 ), 2])
    
40.     }
    
41.     iterN<- iterN+1
    
42.   }
    
43.   
    
44.   
    
45.   #Output of the CMLE.beta funtion
    
46.   Results<- list(beta, iterN, it.log)
    
47.   names(Results) <- c("point estimator", "iteration Num",  "iteration log")
    
48.   return(Results)
    
49. }

复制代码

1. #单个theta的条件极大似然估计；作答矩阵需事先排除满分和零分的被试数据，以及全对和全错的项目数据
   
2. CMLE.theta.1P<-function(resp.vect,D=1.7,beta,theta.start=0,n1=50,tol1=0.0005,...){
   
3.   theta<- theta.start
   
4.   d1.logLikDiff<- 999999999
   
5.   it.log<- matrix(ncol=3,nrow=0)
   
6.   colnames(it.log)<- c("theta_estimator","d1.logLik","d2.logLik")
   
7.   iterN=1 #iteration number
   
8.   
   
9.   while(d1.logLikDiff>tol1 & iterN<n1){
   
10.     probs.1<- 1/(1+exp(-D*(theta-beta)))
    
11.     probs.0<- 1-probs.1
    
12.    
    
13.     #First derivative
    
14.     d1.logLik<- sum(D * (resp.vect - probs.1))
    
15.    
    
16.     #Second derivative
    
17.     d2.logLik<- -D*sum(probs.1 * probs.0)
    
18.     it.log<- rbind(it.log,c(theta,d1.logLik,d2.logLik))
    
19.    
    
20.     #Newton-Raphson
    
21.     theta<- theta-(d1.logLik / d2.logLik)
    
22.    
    
23.     #Convergence critieria
    
24.     if(nrow(it.log) > 1){
    
25.       d1.logLikDiff<- abs(d1.logLik - it.log[(nrow(it.log) - 1), 2])
    
26.     }
    
27.     iterN<- iterN+1
    
28.   }
    
29.   
    
30.   #Standard error
    
31.   SE.theta<- 1/sqrt(-1*d2.logLik)
    
32.   
    
33.   #Output of the CMLE.theta function
    
34.   Results <- list(theta,SE.theta,iterN,it.log)
    
35.   names(Results) <- c("point estimator","standard error", "iteration Num", "iteration log")
    
36.   return(Results)
    
37. }

复制代码

calsunny 发表于 2014-12-2 14:29
Resp.Matrix 和 data1$U_Matrix 是同样的吗？

CMLE.theta.1P（）和CMLE.beta.1P(） 这两个函数不在，所 ...

Resp.Matrix是参数，data1$U_Matrix是具体的作答矩阵。麻烦你帮我看看，谢谢了！

超级大菜鸟 发表于 2014-12-3 12:53
Resp.Matrix是参数，data1$U_Matrix是具体的作答矩阵。麻烦你帮我看看，谢谢了！

OK. I will. But it takes time.

Do you have any particular reasons for not using Winstep or Bilog but writing the functions by yourself?

Is Resp.Matrix 是具体的作答矩阵 while data1$U_Matrix is 参数, right?

Could you double check where data1$U_Matrix is from?

I believe the problem is here:
theta.temp[i]<- cmle.theta.temp[[i]]> JMLE.1P(data1$U_Matrix)

I did not get where data1$U_Matrix is from.

This might be the reason you got the error message:
Error in JMLE.1P(Resp.Matrix = resp.vect, n2 = 20, tol2 = 0.005) :
  object 'theta.temp' not found

calsunny 发表于 2014-12-4 04:27
Is Resp.Matrix 是具体的作答矩阵 while data1$U_Matrix is 参数, right?

Could you double check where ...

data1$U_Matrix是我事先就模拟出来的作答矩阵。
之前忘了发模拟这个矩阵的函数。

1. Resp.Matrix.1P<- function(N,m,D=1.7) {
   
2.   #True value of Difficulty
   
3.   b_T<- rnorm(m)
   
4.   for (j in 1:m) {
   
5.     while (b_T[j] < (-3) | b_T [j] > 3) b_T[j]<- rnorm(1)
   
6.   }
   
7.   
   
8.   #True value of Theta
   
9.   theta_T<- rnorm(N)
   
10.   for (i in 1:N) {
    
11.     while (theta_T[i] < (-3) | theta_T[i] > 3) theta_T[i]<- rnorm(1)
    
12.   }
    
13.   
    
14.   #Rand_Matrix: uniform distribution u~(0,1)
    
15.   #P_Matrix: probability matrix
    
16.   #U_Mattix: generated response matrix
    
17.   Rand_Matrix<- matrix(runif(N*m,0,1),nrow=N,ncol=m)
    
18.   P_Matrix<- matrix(0,N,m)
    
19.   U_Matrix<- matrix(0,N,m)
    
20.   for (i in 1:N) {
    
21.     for(j in 1:m) {
    
22.       P_Matrix[i,j]<- 1/(1 + exp(-D * (theta_T[i] - b_T[j])))
    
23.       if (Rand_Matrix[i,j] <= P_Matrix[i,j]) U_Matrix[i,j]<- 1
    
24.       else U_Matrix[i,j]<- 0  
    
25.     }
    
26.   }  
    
27.   Results<- list(N<-N,m<-m,b_T<-b_T,theta_T<-theta_T,P_Matrix<-P_Matrix,U_Matrix<-U_Matrix)
    
28.   names(Results)<- c("N","m","b_T","theta_T","P_Matrix","U_Matrix")
    
29.   return(Results)
    
30. }

复制代码

data1就是这个函数的结果

1. data1<- Resp.Matrix.1P(4000,100)

复制代码

calsunny 发表于 2014-12-4 03:27
OK. I will. But it takes time.

Do you have any particular reasons for not using Winstep or Bi ...

我是想试着练习一下，方便以后进一步学习IRT。

calsunny 发表于 2014-12-4 03:27
OK. I will. But it takes time.

Do you have any particular reasons for not using Winstep or Bi ...

我是想试着练习一下，方便以后进一步学习IRT。