学习IRT或MATLAB或江西师大的同学们，请来帮个忙！

不知能否给一个利用MATLAB模拟作答矩阵的代码，就是带入irt模型进行MCMC模拟的那种。谢谢啦。

别吓人啊。你知道什么叫MCMC模拟吗？
MC模拟的代码，我们手上还有一些；MCMC可就没谁懂了。

别吓人啊。你知道什么叫MCMC模拟吗？
MC模拟的代码，我们手上还有一些；MCMC可就没谁懂了。

笔为剑 发表于 14-4-11 16:12
别吓人啊。你知道什么叫MCMC模拟吗？
MC模拟的代码，我们手上还有一些；MCMC可就没谁懂了。

马尔科夫蒙特卡洛.....

模拟作答矩阵哪里需要用到马尔科夫链？只要很简单地套作答概率公式就行。
代码可以给你，就看你会不会用了。

function [Sita,A,B,U]=MonteCarloSimulation( )

% N is the number of examinees
% m is the number of items
%
% Sita is the array of the ability parameter,size:1*N.
% A is the array of the discrimination parameter,size:1*m.
% B is the array of the difficulty parameter,size:1*m.
N=200;
m=100;

% 先模拟能力参数,先验分布为标准正态分布
rand('seed',sum(100*clock));
randn('seed',sum(100*clock));
for i=1:N
    flag=1;              % 置标志位为1
    while (flag==1)
        t1=randn;        % 调用库函数randn,用于生成正态分布随机数
        if ((t1<=3)&(t1>=-3))    % 控制范围
            flag=0;      % 得到能力范围内的值,就将flag置为0
        end
    end
    Sita(i)=t1;          % t1为中间值
end
% 接着模拟项目参数
% 区分度参数,先验分布为对数标准正态分布
for i=1:m
    flag=1;
    while (flag==1)
        t2=exp(randn);
        if ((t2>=0.2)&(t2<=2.5))  % 控制范围
            flag=0;
        end
    end
    A(i)=t2;
end
% 难度参数,先验分布为标准正态分布
for i=1:m
    flag=1;
    while (flag==1)
        t3=randn;
        if ((t3<=3)&(t3>=-3))     % 控制范围
            flag=0;
        end
    end
    B(i)=t3;
end
%模拟产生得分阵
for i=1:N-1  
    flag=0;
    while(flag==m|flag==0)  %当生成的被试反应异常（都没答对或都答对）时，则重生成
        flag=0;
        for j=1:m   
            P(i,j)=1/(1+exp(-1.702*A(j)*(Sita(i)-B(j))));   % 根据2PLM计算能力为Sita(i)的被试答对项目j的概率
            r=rand;                     % 生成一[0,1]之间的随机小数
            if P(i,j)>=r
                U(i,j)=1;               % 生成得分矩阵
                flag=flag+1;  
            else U(i,j)=0;
            end
        end
    end
end
flag1=sum(U);
flag=0;
%下面产生第N行
for k=1:m
    if(flag1(k)==0|flag1(k)==N-1) %对异常项目（没有人答对或所有人都答对的项目）最后一个被试的反应能够取反（）
        if(flag1(k)==0)
            U(N,k)=1;flag=flag+1;
        else
            U(N,k)=0;         
        end
    else              % 正常项目对应最后一个被试的反应随机产生
        P(N,k)=1/(1+exp(-1.702*A(k)*(Sita(N)-B(k))));   % 根据2PLM计算能力为Sita(i)的被试答对项目j的概率
        r=rand;                     % 生成一[0,1]之间的随机小数
        if P(N,j)>=r
            U(N,j)=1;               % 生成得分矩阵
            flag=flag+1;   
        else U(N,j)=0;
        end
    end
end
%保证第N个项目正常
if flag==0               
    U(N,m)=1;
end
if flag==m
    U(N,m)=0;
end
end
  U

                    
   
        

首先谢谢您。基本看懂了，不过有个地方
：“ if P(N,j)>=r
            U(N,j)=1;               % 生成得分矩阵
            flag=flag+1;   
        else U(N,j)=0;”
从一个概率值转化为0-1值，为什么与一个随机数比较，不应该是代入贝努力函数生成吗？虽然结果可能一样，但是道理好像说不通。

杨春鑫 发表于 14-4-11 16:26
首先谢谢您。基本看懂了，不过有个地方
：“ if P(N,j)>=r
            U(N,j)=1;               % 生成 ...

因为被试在作答时候有可能猜测也有可能会失误，带有一定随机性。
蒙特卡洛模拟都是这样模拟作答的。