preserve.
set printback=off.
*-------------------------------------------------------------.
* 数量化３類マクロ for spss. by k.hori(hori@ec.kagawa-u.ac.)
*    (欠損値のない場合）                      1996-03-18.
*-------------------------------------------------------------.
define hayasi3 (ivar=!charend('/') 
/maxis=!default(5)!charend('/')
/pscore=!default(NO)!charend('/')
/ascore=!default(NO)!charend('/')
/missing=!default(NO)!charend('/')
/sdtype=!default(3)!charend('/')
/each=!default(0)!charend('/')
).
preserve.
set printback=off.
*---------------------------------------------------------------------------
* 欠損値データの処理.
*---------------------------------------------------------------------------.
*元のデータを保存.
save outfile='h3_tmp1.sav'.
!let !miss2=!upcase(!missing).
!if (!miss2='') !then .
!let !miss2='NO'.
!ifend.
!if (!miss2 <> 'NO') !then.
count ms_=!ivar (missing).
select if ms_=0.
!ifend.
!let !nvs1=''.
!let !nvar=!length(!nvs1).
set mxloops=40000.
!let !ascore2=!upcase(!ascore).
!if (!ascore2='') !then .
!let !ascore2='NO'.
!ifend.
*****************************************************************.
matrix.
compute colname={" ","******    ","数量化３","類","******"," ",
"v.0.20 "}.
print colname/format=a10.
compute nprint=8.
*design の出力を1,2,3...とするため num.
compute num={1,2,3,4,5,6,7,8,9,10,11,12,13,14,15}.
compute num=t(num).
compute num2={'1','2','3','4','5','6','7','8','9','10','11','12','13','14','15'}.
compute nx={'error'}.
*get x /file=* /var=!ivar /missing=0/names=nx.
get x /file=* /var=!ivar /sysmis=omit/names=nx.
compute nsample=nrow(x).
* print nsample.
!if (nsample !gt 0) !then. 
*変数数 #nvar.
compute #nvar=ncol(nx).
!let !varname=''.
* cmaxs 各itemの最大カテゴリー数.
compute cmaxs=cmax(x).

print nsample.
*カテゴリの１，０展開 結果を z に.
compute z=x(:,1).
loop i=1 to #nvar.
compute x1=x(:,i).
compute x1={num(1:cmaxs(1,i),1);x1}.
compute d=design(x1).
compute #row1=nsample+cmaxs(i).
compute d=d((cmaxs(i)+1):#row1,1:cmaxs(i)).
compute z={z,d}.
end loop.
*カテゴリ総数 ctotal.
compute ctotal=ncol(z)-1.
compute x=z(:,2:(ctotal+1)).
*各カテゴリの反応数.
compute ncat=csum(x).
compute tres=msum(x).
compute ntotal=msum(ncat).
release z,x1,d.
*---------------------------------------------------.
* Compute the pooled within groups SSCP matrix.
* Compute the beteen gropus SSCP matrix.
*---------------------------------------------------.
compute sumr=sqrt(rsum(x)).
compute x2=x.
loop i=1 to ctotal.
compute x2(:,i)=x2(:,i)/sumr.
end loop.
compute spt=sscp(x2).
compute sumr=rsum(x).
release x2.
*--------------------------------.
* 特異値分解.
*--------------------------------.
compute d=mdiag(ncat).
compute dr_inv=inv(sqrt(d)).
compute w=dr_inv*spt*dr_inv.
compute label={'total',num2}.
call svd(w,evc,evl,v).
compute evc=dr_inv*evc.
compute evl=diag(evl).
compute evl=evl(2:nrow(evl),1).
*----------------------------------.
*各種寄与率.
*----------------------------------.
compute maxis=!maxis.
compute j=1.
compute evl2=evl.
compute evl3=evl.
compute q=#nvar/(#nvar-1).
print t(evl)/title "固有値".
loop  if (evl(j,1)>(1/#nvar)).
compute evl2(j,1)=(q*(evl(j,1)-1/#nvar))**2.
compute evl3(j,1)=(q*(sqrt(evl(j,1))-1/#nvar))**2.
compute j=j+1.
end loop.
print (1/#nvar)/title="許容最小固有値（１／アイテム数)(Nishisato,1994)"
  /format=f9.7.

do if (maxis >= j) .
compute maxis=j-1.
end if.
compute j=j-1.
compute evl2sum=msum(evl2(1:j,1)).
compute evl3sum=msum(evl3(1:j,1)).
compute j=maxis.
compute evl2=evl2(1:maxis,1).
compute evl3=evl3(1:maxis,1).
compute evl4=evl(1:maxis,1)&/sqrt(mssq(evl)).
compute evl5=evl.
*.
loop j1=1 to nrow(evl).
   compute kk=0.
 loop j2=1 to j1.
   compute kk=kk+evl(j2,1)*evl(j2,1).
 end loop.
   compute evl5(j1,1)={sqrt(kk)&/sqrt(mssq(evl))}.
   compute j3=evl5(j1,1).

   do if (j1>1).
     compute evl5(j1,1)={j3-j4}.
   end if.
   compute j4=j3.
end loop.
*.
compute cnam={"固有値","ratio","Benzecri","  ratio","Greenacr","e  ratio",
    "岩坪","岩坪２"}.
compute xx={evl(1:maxis,1),evl(1:maxis,1)&/msum(evl),evl2,(evl2&/evl2sum),
            evl3(1:maxis,1),(evl3(1:maxis,1)/evl3sum),
            evl4(1:maxis,1),(evl5(1:maxis,1)&/msum(evl5))}.
print xx/cnames=cnam/rnames=num2/title="寄与率"/format=f8.3.
release w,xx.
*reliablity.
compute alpha=make(maxis,1,#nvar).
compute alpha=(alpha-make(maxis,1,1)&/(evl(1:maxis,1)))&/(#nvar-1).
*chi-square.
compute chi=make(maxis,1,((nsample+ctotal-#nvar-1)/2-ntotal+1)).
compute chi=chi&*ln(make(maxis,1,1)-(evl(1:maxis,1))).
compute df=make(maxis,1,nsample+ctotal-#nvar).
compute prob=make(maxis,1,0).
loop i=1 to maxis.
compute df(i,1)=df(i,1)-2*i-1.
compute prob(i,1)=1-chicdf(chi(i,1),df(i,1)).
end loop.
print {alpha,chi,df,prob}/title="軸の統計量"
      /format=f12.4/rnames=num2/clabels="信頼性α","χ2","df","p".
compute evc=evc(:,2:(maxis+1)).
*------------------.
*カテゴリーラベル.
*------------------.
compute catlab={nx(1,1),num2(1,2:cmaxs(1,1))}.
do if #nvar > 1.
loop i=2 to #nvar.
compute catlab={catlab,nx(1,i),num2(1,2:cmaxs(1,i))}.
end loop.
end if.
*------------------------------------.
*固有ベクトルの調整＝＞数量化の係数.
*------------------------------------.
compute dcoef=evc.
*print dcoef/format=f7.2/title="eigen vector"
* /rnames=catlab/cnames=num2.
*平均値・標準偏差.
!if (!sdtype = '3') !then.
compute sd=make(1,maxis,1/sqrt(ntotal)).
compute sd=sd&*t(sqrt(evl(1:maxis,1))).
compute means=make(1,maxis,0).
!ifend.

*平均・標準偏差終了------------------.

loop i=1 to ctotal.
compute dcoef(i,:)=(dcoef(i,:)-means)&/sd.
end loop.
compute clab={"freq",num2}.
*print {t(ncat),dcoef}/format=f7.2/title="駒澤型数量化３類　カテゴリ数量"
*/rnames=catlab/cnames=clab.
compute coef=dcoef.
compute  evalue=mdiag(t(evl(1:maxis,1))).
compute coef=coef*evalue.
print {t(ncat),coef}/format=f7.2/
      title="HOMALS型 対応分析型　カテゴリ数量（正準正規化・主軸型・主座標）"
      /rnames=catlab/cnames=clab.
compute eta2=t(t((coef)&**2)*d)&/nsample.
compute inertia=diag(inv(d)*spt).
compute dcoef=t(dcoef).
compute coef=t(dcoef).
compute coef={t(ncat),coef}.
compute clab={"freq",num2}.
*print coef/format=f6.2/title="　　　　カテゴリ数量"
*  /rnames=catlab/cnames=clab.
*================.
*  各変数の得点.
*================.
compute loading1=make(ctotal,1,nsample).
compute loading1=sqrt((loading1-t(ncat))&/t(ncat))).
compute loading=loading1.
compute k=1.
loop if (k< maxis).
  compute loading={loading,loading1}.
  compute k=k+1.
end loop.
compute loading=t(dcoef)&/loading.
compute loading=loading*evalue.
print loading/format=f6.2/
      title="負荷量"
      /cnames=num2/rnames=catlab.
compute loading1=loading&**2.
print loading1/format=f6.2/
      title="相関（負荷量）の２乗（対応分析）"
      /cnames=num2/rnames=catlab.
print eta2/format=f6.3/
      title="軸×カテゴリーごとのinertia(分散)"
      /cnames=num2/rnames=catlab.
compute ainertia=eta2*inv(evalue)&/#nvar.
print ainertia/format=f6.3/
      title="inertia(固有値)への絶対的寄与率"
      /cnames=num2/rnames=catlab.
compute tinertia=rsum(eta2)&/rsum(loading1).
compute summary={rsum(loading1),(rsum(eta2)&/msum(eta2)),
       (tinertia&/msum(tinertia))}.
print summary 
    /format=f9.3/ title="集約統計 (対応分析型)"
      /clabels="Quality","inertia","全inertia"/rnames=catlab.
*********************************.
compute varcorr=make(#nvar,maxis,0).
compute range=make(#nvar,maxis,0).
compute sk=0.
compute rowname={num2}.
loop i=1 to #nvar.
compute ek=sk+cmaxs(1,i).
compute sk=sk+1.
compute inertia=eta2(sk:ek,:).
compute varcorr(i,:)=(csum(inertia)).
compute m=t(dcoef(:,sk:ek)).
compute range(i,:)=(cmax(m)-cmin(m)).
compute sk=ek.
end loop.
print varcorr/format=f6.2/
      title="相関係数の２乗（HOMALS: discriminant measure）"
      /cnames=num2/rnames=nx.
*print range/format=f6.2/
*      title="範囲 （数量化３類: range:駒澤型）"
*      /cnames=num2/rnames=nx.
*print (range*sqrt(evalue))/format=f6.2/
*      title="範囲 （数量化３類: range:正規化）"
*      /cnames=num2/rnames=nx.
print (range*evalue)/format=f6.2/
      title="範囲 （数量化３類: range:正準正規化）"
      /cnames=num2/rnames=nx.
*============.
* 個体数量.
*============.
!let !pscore2=!upcase(!pscore).
!if (!pscore2='') !then.
!let !pscore2='NO'.
!ifend.
!if ((!pscore2 <> 'NO') !or (!ascore2 <> 'NO')) !then.
compute dscore=x*t(dcoef).
loop i=1 to maxis.
compute dscore(:,i)=dscore(:,i)&/sumr.
end loop. 
compute label={num2(1,1:maxis)}.
!ifend.
!if (!pscore2 <> 'NO') !then .
print dscore/title="個体数量（標準得点）"/cname=label/format=f8.4.
!ifend.
*------------------------------------.
* 個体数量のカレントファイルへの出力.
*------------------------------------.
!if (!ascore2 <> 'NO') !then.
compute vlabel={"h3s1","h3s2","h3s3","h3s4","h3s5","h3s6","h3s7","h3s8","h3s9","h3s10","h3s11","h3s12","h3s13","h3s14","h3s15","h3s16"}.
save dscore/outfile="h3_score.sav"/names=vlabel.
!ifend.
end matrix.
!if (!ascore2 <> 'NO') !then .
match files file="h3_score.sav"/file=*/map.
exec.
!ifend.
!else .
end matrix.
!ifend.
restore.
!enddefine.
restore.

/*hayasi3 ivar=A1 a2 b3 /pscore=yes/missing=NO.*/
/*show  workspace*/.
/*set workspace 14000*/.
/*hayasi3 ivar=t30001 to t30021/ascore=yes /missing=yes*/.