包含在主程式 main.m 的自建函數檔
一、波氏演算法(BaumWelch_Algorithm)--參數估算 function
[alpha,beta,app_A,app_B,app_pi,prob1]=BaumWelch_Algorithm(O,A,B,pi)
% (ㄧ) 算出 alpha,beta,P(Ow|λ)值
%% 正算程序
[alpha,probf]=Forward_procedure(O,A,B,pi);
%% 逆向算程序
[beta,probb]=Backward_procedure(O,A,B,pi);
%% 正算和逆算程序
[probfb]=FBward_procedure(O,A,alpha,beta);
% (二) GHMM 參數重估演算法,先定義兩個參數 phi,rho [phi,rho]=arg2(O,A,B,alpha,beta,probfb);
%% 3.符號狀態機率變動法
[app_A,app_B,app_pi]=em3(O,B,phi,rho);
% 由(二)GHMM 參數重估演算法之 3.符號狀態機率變動法,可產生我們所估來的隨 時間變動的 A 矩陣,B 矩陣和 pi 向量,可算出新的機率值((probw=連乘 w 的 P(Ow|
λnew))),取 rhog,用 rhog(sumw(P(Ow|λnew)))
[alpha,probf]=Forward_procedure(O,app_A,app_B,app_pi);
prob1=sum(log(probf));
二、 正算程序--算出 P(Ow|λ)值
function [alpha,probf]=Forward_procedure(O,A,B,pi)
% 1.定初始值 for w=1:size(O,1) for i=1:3
alpha(w,i,1)=pi(i)*B(i,O(w,1),1,w);
end end
% 2.產生遞迴正算變數 for t=1:size(O,2)-1 for j=1:3
for w=1:size(O,1) sumaA=0;
for i=1:3
sumaA=sumaA+alpha(w,i,t)*A(i,j,t);
end
alpha(w,j,t+1)=sumaA*B(j,O(w,t+1),t+1,w);
end end end
% 3.終止條件--probf=P(Ow|λ) for w=1:size(O,1)
probf(w)=0;
for i=1:3
probf(w)=probf(w)+alpha(w,i,size(O,2));
end
三、逆算程序--算出 P(Ow|λ)值
function [beta,probb]=Backward_procedure(O,A,B,pi)
% 1.定初始值 for i=1:3
for w=1:size(O,1)
beta(w,i,size(O,2))=1;
end end
% 2.產生遞迴逆算變數 for t=size(O,2)-1:-1:1 for i=1:3
for w=1:size(O,1) sumABb=0;
for j=1:3
sumABb=sumABb+A(i,j,t)*B(j,O(w,t+1),t+1,w)*beta(w,j,t+1);
end
beta(w,i,t)=sumABb;
end end end
% 3.終止條件--probb=P(Ow|λ) for w=1:size(O,1)
sumpBb=0;
for i=1:3
sumpBb=sumpBb+pi(i)*B(i,O(w,1),w)*beta(w,i,1);
end
probb(w)=sumpBb;
end
四、正算和逆算程序--算出 P(Ow|λ)值
function [probfb]=FBward_procedure(O,A,alpha,beta ) [W,T]=size(O);
for w=1:W for t=1:T sumab=0;
for i=1:3
sumab=sumab+alpha(w,i,t)*beta(w,i,t);
end end
probfb(w)=sumab;
end
五、先定義兩個參數 fa(w,i,j,t)=P(s(w,t)=i,s(w,t+1)=j|Ow, λ);lo(w,i,t)=P(s(w,t)=i|Ow,λ)
function [phi,rho]=arg2(O,A,B,alpha,beta,probf) [W,T]=size(O);
for w=1:W for i=1:3
for t=1:T-1
phi(w,i,j,t)=alpha(w,i,t)*A(i,j,t)*B(j,O(w,t+1),t+1,w)*beta(w,j,t+1)/
probf(w);
end end end end
for w=1:W for i=1:3 for t=1:T
rho(w,i,t)=alpha(w,i,t)*beta(w,i,t)/probf(w);
end end end
六、符號狀態機率變動法
% 估出 A,B,pi 的估計值 app_A,app_B,app_pi,且可檢驗矩陣的列和等於一,值的 小數點位數取到小數後第三位,且去除等於零的機率值(ps:app_A 取三列的平均 機率;且 0<pi2=a12<0.25,0<pi3=a13<0.05)
function [app_A,app_B,app_pi]=em22(O,B,phi,rho) [W,T]=size(O);
app_A=zeros(3,3,T-1);app_B=zeros(3,3,T,W);app_pi=zeros(1,3);
for t=1:T-1 for j=1:3 sumaphi=0;
sumarho=0;
for w=1:W for i=1:3
sumaphi=sumaphi+phi(w,i,j,t);
sumarho=sumarho+rho(w,i,t);
end end
if sumarho==0 sumarho=0.001;
end
for i=1:3
app_A(i,j,t)=sumaphi/sumarho;
end end end
% 每題組(兩題為一題組)間因相互獨立，因此題組間 app_A 均相等,取 app_A 三 列的平均機率
app_A1=app_A;
for t=2:2:T-1
app_A(:,:,t)=mean(app_A1(:,:,:),3);
end
%%% 可檢驗 app_A 矩陣的列和等於一,即 sum(app_A,2)=1 for i=1:3
for t=1:T-1
app_A(i,3,t)=0.001;
app_A(i,1,t)=app_A(i,1,t)-0.001;
end end end
for t=1:T-1 for i=1:3
if app_A(i,2,t)<app_A(i,3,t) a12=app_A(i,2,t);
a13=app_A(i,3,t);
app_A(i,2,t)=a13;
app_A(i,3,t)=a12;
end end
end for t=1:T for m=1:3 for j=1:3
sumbrho22=0;
for w=1:W if m==1
[ww,tt]=find(O==2 | O==3);
rhorho=rho;
for k=1:size(ww,1)
rhorho(ww(k),j,tt(k))=0;
end
newrho=sum(rhorho,1);
end if m==2
[ww,tt]=find(O==1 | O==3);
rhorho=rho;
for k=1:size(ww,1)
rhorho(ww(k),j,tt(k))=0;
end
newrho=sum(rhorho,1);
end if m==3
[ww,tt]=find(O==1 | O==2);
rhorho=rho;
for k=1:size(ww,1)
rhorho(ww(k),j,tt(k))=0;
end
newrho=sum(rhorho,1);
end
sumbrho22=sumbrho22+rho(w,j,t);
end
if sumbrho22==0 sumbrho22=0.001;
end
app_B(j,m,t,w)=newrho(1,j,t)/sumbrho22;
end end end
end
%%% 可檢驗 app_B 矩陣的列和等於一,即 sum(app_B,2)=1 for w=1:W
for t=1:T
app_B(:,3,t,w)=[0 0 1];
app_B(:,2,t,w)=1-app_B(:,1,t,w)-app_B(:,3,t,w);
end end
% 要使 app_B 矩陣裡的猜測答對符號(即為 c)<0.25 for w=1:W
for t=1:T
if app_B(2,1,t,w)>0.25
kk=0.1875+0.0625*rand(1);
app_B(2,2,t,w)=app_B(2,2,t,w)+(app_B(2,1,t,w)-kk);
app_B(2,1,t,w)=kk;
end end end
for i=1:3 sumrho=0;
for w=1:W
sumrho=sumrho+rho(w,i,1);
end
app_pi(i)=sumrho/W;
end
%%% 可檢驗 app_pi 矩陣的列和等於一,即 sum(app_pi,2)=1 if app_pi(3)==0
app_pi(3)=0.001;
app_pi(1)=app_pi(1)-0.001;
end
if app_pi(2)<app_pi(3) pi2=app_pi(2);
pi3=app_pi(3);
app_pi(2)=pi3;
app_pi(3)=pi2;
end
七、KN-IRT 估計法
function [KN_Ptheta,X,X1,est,KN_theta1,n]=KNIRT1(O) [W,T]=size(O);
X=zeros(W,T+1); % X 表受試者之測驗得分情形 rt=zeros(T+1); % rt 表試題 t 之「相關鑑別指數」
U=zeros(W,1); % U 表相關加權測驗總分 KN_Ptheta=zeros(W,T); %預設答對機率值矩陣 for w=1:W
if O(w,t)==1 % O 表受試 w 作答試題 t 答對與否之指示 變數
X(w,t)=1;
X(w,T+1)=X(w,T+1)+1;
else
X(w,t)=0;
end end
end for t=1:T
rt(t)=corr2(X(:,t),X(:,T+1));
u(t)=(1+rt(t))/2;
end
for w=1:W for t=1:T
U(w)=U(w)+u(t)*X(w,t); % U 表相關加權測驗總分 end
end
KN_theta=zscore(U); %受試者 w 之能力參數估計值 [m,n]=sort(KN_theta);
est_theta=zeros(W,1);
for w=1:W
X1(w,:)=X(n(w),:);
KN_theta1(w,:)=KN_theta(n(w),:);
end
est=norminv(1/(W+1):1/(W+1):W/(W+1),0,1);
%計算受試能力為 real_theta 者之試題 t 之「核平滑化無參數試題特徵曲線」
for t=1:T
for j=1:length(est) dem=0;
num=0;
for w=1:W
dem=dem+exp(-(W^0.4*(KN_theta(w)-est(j)).^2)/2.42);
num=num+exp(-(W^0.4*(KN_theta1(w)-est(j)).^2)/2.42).*X1(w,t);
end
KN_Ptheta(j,t)=num/dem; %答對機率函數 end
end
八、GHMM 與 KN-IRT 之整合估計 function
[GKN_Ptheta,GKN_theta1,gkn_t,gkn_p]=GHMM_KNIRT1(pi,A,c,KN_Ptheta,X1,e st)
[W,T]=size(KN_Ptheta);
U=zeros(W,1);
GKN_Ptheta=zeros(W,T); %已分離猜測及未答之無參數 IRT 特徵曲線
%計算已分離猜測及未答之無參數 IRT 特徵曲線 for w=1:W
GKN_Ptheta(w,1)=(KN_Ptheta(w,1)-pi(2)*c(1))/pi(1);
for t=2:T for w=1:W
GKN_Ptheta(w,t)=(KN_Ptheta(w,t)-A(1,2,t-1)*c(t))/A(1,1,t-1);
end end
%求已分離猜測及未答之無參數 IRT 受試者 w 之能力參數估計值 for w=1:W
U(w)=pi(1)*GKN_Ptheta(w,1);
for t=2:T
U(w)=U(w)+A(1,1,t-1)*GKN_Ptheta(w,t);
end end
GKN_theta1=zscore(U);
%已分離猜測及未答之無參數 IRT 特徵曲線遞迴重估
S=1;gk=zeros(100,1);count=1; %count 表遞迴重估次數 while S > 0.001
GKN_Ptheta=zeros(W,T); %已分離猜測及未答之無參數 IRT 特徵 曲線
for t=1:T
for w=1:length(est) dem=0;
num=0;
for j=1:W
num=num+exp(-(W^0.4*(GKN_theta1(j)-est(w)).^2)/2.42).*X1(j,t);
dem=dem+exp(-(W^0.4*(GKN_theta1(j)-est(w)).^2)/2.42);
end
GKN_Ptheta(w,t)=num/dem;
end end for w=1:W
GKN_Ptheta(w,1)=(GKN_Ptheta(w,1)-pi(2)*c(1))/pi(1);
end
for t=2:T for w=1:W
GKN_Ptheta(w,t)=(GKN_Ptheta(w,t)-A(1,2,t-1)*c(t))/A(1,1,t-1);
end end
for t=1:T for w=1:W
if GKN_Ptheta(w,t) >= 1 GKN_Ptheta(w,t)=0.9999;
end
if GKN_Ptheta(w,t) <= 0 GKN_Ptheta(w,t)=0.0001;
end end end
gkn_p(:,:,count)=GKN_Ptheta;%%測試用 U=zeros(W,1);
U(w)=pi(1)*GKN_Ptheta(w,1);
for t=2:T
U(w)=U(w)+A(1,1,t-1)*GKN_Ptheta(w,t);
end end
GKN_theta2=zscore(U); %下一次遞迴重估 theta V=0;
gkn_t(:,:,count)=GKN_theta1;%%測試用 for w=1:W
V=V+(GKN_theta2(w)-GKN_theta1(w))^2;
end S=V/W;
GKN_theta1=GKN_theta2;
gk(count)=S;
count=count+1;
end
%計算已分離猜測及未達之無參數 IRT 特徵曲線 P(theta) for t=1:T
for w=1:W dem=0;
num=0;
for j=1:W
num=num+exp(-(W^0.4*(GKN_theta1(j)-GKN_theta1(w)).^2)/2.42).*X1(j,t);
dem=dem+exp(-(W^0.4*(GKN_theta1(j)-GKN_theta1(w)).^2)/2.42);
end
GKN_Ptheta(w,t)=num/dem;
end end for w=1:W
GKN_Ptheta(w,1)=(GKN_Ptheta(w,1)-pi(2)*c(1))/pi(1);
end
for t=2:T for w=1:W
GKN_Ptheta(w,t)=(GKN_Ptheta(w,t)-A(1,2,t-1)*c(t))/A(1,1,t-1);
end end for t=1:T for w=1:W
if GKN_Ptheta(w,t) >= 1 GKN_Ptheta(w,t)=0.9999;
end
if GKN_Ptheta(w,t) <= 0 GKN_Ptheta(w,t)=0.0001;
end end end
九、估計作答策略變數矩陣
%% 3.估出新的 y:當 y0(w,t,3)=0,則 y(w,t,1)+y(w,t,2)=1 估出新的 y(w,t,1)
function y=e_y(O,U,y0,pi,c,Ptheta) [W,T]=size(O);
for w=1:W for t=1:T
if y0(w,t,3)==0 y(w,t,3)=0;
y(w,t,1)=(pi(1,1)*Ptheta(w,t)^U(w,t)*(1-Ptheta(w,t))^(1-U(w,t)))/(pi(
1,1)*Ptheta(w,t)^U(w,t)*(1-Ptheta(w,t))^(1-U(w,t))+pi(1,2)*c(t)^U(w,t )*(1-c(t))^(1-U(w,t)));
y(w,t,2)=(pi(1,2)*c(t)^U(w,t)*(1-c(t))^(1-U(w,t)))/(pi(1,1)*Ptheta(w, t)^U(w,t)*(1-Ptheta(w,t))^(1-U(w,t))+pi(1,2)*c(t)^U(w,t)*(1-c(t))^(1-U(w,t)));
else
y(w,t,3)=1;
y(w,t,1)=0;
y(w,t,2)=0;
end end end
%%% 可檢驗 當 y(w,t,3)=0,則 y(w,t,1)+y(w,t,2)=1;sum(y,3)=1