Attribute VB_Name = "polycorr11" '使用法: '(a)相関係数だけを求めるときは関数 polycorr を使うと求めることができる。=polycorr(範囲指定) '(b)いろいろな指標も求める場合はマクロ, poly_corr を起動して範囲を指定する。 '(c)optionを一時的に変えたいときはマクロ,poly_corr_option を起動してオプションを指定し,範囲を指定する。 '(d)option を変えてしまいたいときは sub init0 のなかの変数の値を書き換える。 '範囲はいずれもクロス表を指定する。マウスを使うこともできる。 'poly_corr_Matrix_from_raw_data は素データ行列から多分相関係数行列を求める。 '(a)データ以外のセルからマクロ poly_corr_Matrix_from_raw_data を起動する。 '(b)データは変数名が最初についていればそれを使うなければ,勝手に変数名をつける。 '(c)データ範囲の指定はデータの下および右がempty ならば左上のセルを指定すればいい。 '(d)データの下がempty なら,一列(変数名またはデータ)を指定すればいい。 '(e)すべてのデータの範囲を指定してもいい。 '(f)データに欠損値がないようにあらかじめ処理しておくこと。 'vba に移植してからいくつかチェックをいれた。 '(1)2つの推定法のML(最尤)法と最小χ2法ではχ2法の推定値に納得いかない推定値があり,ML法を使うのがいいであろう。 '(2)2-step 推定法とjoint 推定法ではjoint 推定法は収束がはっきりと遅い場合が多くある。10倍以上時間がかかったりする。2-step推定法がいいであろう。 '(3)cdf のルーティンとして,Uebersax氏が最初から入れているのは3つであるが,1つは落とした。alnorm, NPROBをそのまま入れている。ほかにexcel のnormsdist を追加した。excel のnormsdist は最速というわけでもない。alnormは精度に少し問題がある。NPROBを使うのがいいであろう。 '(4)χ2値のp値を求めるルーティンはexcel のchidist に差し替えている。 '(5)メタセルは使用できない。 DefDbl A-H, O-Z DefInt I-N 'common 変数 Private icode Public DataRange As Range '/FASTER/ Private JVARY '/ICALL/ Private ICALLS Private IMAX Private nconv 'STPIT / Private CHISQ Private x(41) Private XMAX(41) Private XMIN(41) Private DELTAX(41) Private ERR(41, 41) Private DELMIN(41) Private MASK(41) Private NV, NTRACE, MATRIX '/BLOCK1/ Private TITLE As String Private VARLBL(8, 2) Public IOP(15) Private FO(364) Private FE(364) Private IVEC(364, 2) Private go1(20) Private GE(20) Private B(81) Private r(81) Private A(1) Private t(2, 18) Private BIAS(2) Private P(2, 81, 18) Private CHISQL, CHISQP, LC, MC, NCELLS, NDIS, NEST, NR, nt, total Private FMARG(3, 18) Private CFMARG(3, 18) Private EMARG(3, 18) Private NEP, SERHO Private SET0(2, 18) Private chisql0, chisqp0, ndf0, corr, nflag, tlimit, nmeta, mc1, mc2, prcorr Private mask2(364) '/seflag/ Private iflagse Const maxindex As Integer = 15 '100 Sub poly_corr_Matrix_from_raw_data() Dim dm As Range Dim row As Integer, col As Integer Dim nc As Integer, x() As Double, y() As Double Dim nx As Integer, k As Integer, j As Integer, kugiri Dim i As Integer Dim xname As String, yname As String Dim n As Integer Dim ccell As Range Dim startcell As String, migicell As String, sitacell As String, endcell As String Dim presult As Variant Set ccell = ActiveCell Set dm = Application.InputBox(prompt:="素データのセル範囲,またはデータの左上端を指定(変数名も可):", TITLE:="順序データの相関表作成", Type:=8) nc = dm.Rows.Count If nc <= 0 Or dm.Columns.Count <= 0 Then MsgBox "範囲または左端セルを指定してください" Exit Sub End If If nc = 1 Then startcell = dm.Address ' MsgBox "s:" + startcell kugiri = InStr(startcell, ":") If kugiri > 0 Then startcell = Left$(startcell, kugiri - 1) migicell = Mid$(dm.Address, kugiri + 1) Else i = 0 Do i = i + 1 Loop Until (Range(startcell).Offset(0, i).Value = Empty) migicell = Range(startcell).Offset(0, i - 1).Address End If i = 0 Do i = i + 1 Loop Until (Range(startcell).Offset(i, 0).Value = Empty) sitacell = Range(startcell).Offset(i - 1, 0).Address endcell = Left$(migicell, InStr(2, migicell, "$") - 1) + Mid$(sitacell, InStr(2, sitacell, "$")) ' MsgBox "e:" + endcell Set dm = Range(startcell + ":" + endcell) End If nc = dm.Rows.Count nx = dm.Columns.Count ' MsgBox nx ' End ccell.Value = "polychoric correlation matrix" Set ccell = ccell.Offset(1, 0) For i = 1 To nx If IsNumeric(dm(1, k)) = False Then ccell.Offset(0, i).Value = dm(1, i) ccell.Offset(i, 0).Value = dm(1, i) Else ccell.Offset(0, i).Value = "V" + Trim(Str(i)) ccell.Offset(i, 0).Value = "V" + Trim(Str(i)) End If Next i ccell.Offset(1, 1).Activate ActiveCell.Value = 1# ActiveCell.Offset(1, 0).Activate ReDim x(nc), y(nc) For k = 2 To nx For j = 1 To k - 1 n = 0 For i = 0 To nc - 1 If i = 0 And IsNumeric(dm(1, k)) = False Then Else x(n) = dm(i + 1, k) y(n) = dm(i + 1, j) n = n + 1 End If Next i presult = csub(x, y, n) ActiveCell.Value = presult If j = k - 1 Then ActiveCell.Offset(0, 1).Activate ActiveCell.Value = 1# ActiveCell.Offset(1, 1 - k).Activate Else ActiveCell.Offset(0, 1).Activate End If Next j Next k ActiveCell.Offset(0, -1).Activate ActiveCell.Value = "N=" ActiveCell.Offset(0, 1).Activate ActiveCell.Value = n MsgBox ("終了") End Sub Function mkindex(x() As Double, ix() As Double, n As Integer) As Integer Dim i As Integer, j As Integer Dim val As Double mkindex = 0 For i = 0 To n - 1 val = x(i) ix(mkindex) = val For j = 0 To mkindex If ix(j) = val Then Exit For End If Next j If j = mkindex Then mkindex = mkindex + 1 End If Next i End Function Sub sortindex(ix() As Double, nx As Integer) Dim i As Integer, j As Integer, im As Integer Dim x As Double, mx As Double For i = 0 To nx - 2 im = i mx = ix(i) For j = i + 1 To nx - 1 If ix(j) < mx Then im = j mx = ix(j) End If Next j If im <> i Then x = ix(i) ix(i) = mx ix(im) = x End If Next i End Sub Function getindex(x As Double, ix() As Double, nx As Integer) As Integer Dim i As Integer For i = 0 To nx - 1 If ix(i) = x Then getindex = i + 1 '**** Exit Function End If Next i End Function Function csub(x() As Double, y() As Double, n As Integer) As Variant Dim t(maxindex, maxindex) As Integer Dim ix(maxindex) As Double, iy(maxindex) As Double Dim nx As Integer, ny As Integer Dim i As Integer, j As Integer, k As Integer csub = "error" nx = mkindex(x, ix, n) ny = mkindex(y, iy, n) If nx = 0 Or ny = 0 Then MsgBox "相異なる変数値が" + CStr(maxindex) + "種類以上ありますので,集計できません" Exit Function End If sortindex ix, nx sortindex iy, ny For i = 0 To n - 1 j = getindex(x(i), ix, nx) k = getindex(y(i), iy, ny) t(j, k) = t(j, k) + 1 Next i csub = polycorr2(t(), nx, ny) End Function Function polycorr2(xx() As Integer, nx, ny) As Variant Dim n As Integer, jj As Integer, i As Integer, j As Integer, ifault As Integer mc1 = nx mc2 = ny If mc1 <= 0 Then polycorr2 = "cell数が0以下です" Exit Function ElseIf mc1 < 2 Or mc2 < 2 Then polycorr2 = "cell数が少ない" Exit Function End If jj = 0 For i = 1 To mc1 For j = 1 To mc2 jj = jj + 1 FO(jj) = xx(i, j) Next j Next i Call init0 Call main ' If nflag = 0 Then polycorr2 = A(1) Else polycorr2 = -A(1) End If If ifault > 0 Then MsgBox (A(1)) If (ifault = 1) Then polycorr2 = "<> Routine did not converge. Try increasing the maximum number of iterations (NITER) in Subroutine TETRA (see Brown, 1977)." Exit Function End If If (ifault = 2) Then polycorr2 = "<> A negative frequency entered, or a row or column marginal is 0 (see Brown, 1977)." End If End If End Function Sub poly_corr_option() Call init0 polycorr1.Show ' Dim xx As Range ' Set xx = Application.InputBox(prompt:="2x2以上の行列指定", TITLE:="polychoric correlation", Type:=8) mc1 = DataRange.Rows.Count mc2 = DataRange.Columns.Count If mc1 <= 0 Then MsgBox ("範囲を指定してください") Exit Sub ElseIf mc1 < 2 Or mc2 < 2 Then MsgBox ("2x2 以上の範囲を指定してください") Exit Sub End If jj = 0 For i = 1 To mc1 For j = 1 To mc2 jj = jj + 1 FO(jj) = DataRange(i, j) Next j Next i icode = 0 ' Call init0 Call main Call display End Sub Function polycorr(xx As Range) As Variant Dim n As Integer mc1 = xx.Rows.Count mc2 = xx.Columns.Count If mc1 <= 0 Then polycorr = "範囲を指定してください" Exit Function ElseIf mc1 < 2 Or mc2 < 2 Then polycorr = "2x2 以上の範囲を指定してください" Exit Function End If jj = 0 For i = 1 To mc1 For j = 1 To mc2 jj = jj + 1 FO(jj) = xx(i, j) Next j Next i Call init0 Call main ' If nflag = 0 Then polycorr = A(1) Else polycorr = -A(1) End If If ifault > 0 Then MsgBox (A(1)) If (ifault = 1) Then polycorr = "<> Routine did not converge. Try increasing the maximum number of iterations (NITER) in Subroutine TETRA (see Brown, 1977)." Exit Function End If If (ifault = 2) Then polycorr = "<> A negative frequency entered, or a row or column marginal is 0 (see Brown, 1977)." End If End If End Function Public Sub poly_corr() Dim xx As Range Set xx = Application.InputBox(prompt:="2x2以上の行列指定", TITLE:="polychoric correlation", Type:=8) mc1 = xx.Rows.Count mc2 = xx.Columns.Count If mc1 <= 0 Then MsgBox ("範囲を指定してください") Exit Sub ElseIf mc1 < 2 Or mc2 < 2 Then MsgBox ("2x2 以上の範囲を指定してください") Exit Sub End If jj = 0 For i = 1 To mc1 For j = 1 To mc2 jj = jj + 1 FO(jj) = xx(i, j) Next j Next i icode = 0 Call init0 Call main Call display End Sub Private Sub display() ' If nflag = 0 Then rz = A(1) Else rz = -A(1) End If 'MsgBox = variant(r) ActiveCell.Value = "polychoric corr = " ActiveCell.Offset(0, 1).Activate ActiveCell.Value = rz ActiveCell.Offset(1, -1).Activate i = NEP + 1 SERHO = Sgn(ERR(i, i)) * (Sqr(Abs(ERR(i, i)))) ActiveCell.Value = "Std err =" ActiveCell.Offset(0, 1).Activate ActiveCell.Value = SERHO ActiveCell.Offset(1, -1).Activate ActiveCell.Value = "Peason corr =" ActiveCell.Offset(0, 1).Activate ActiveCell.Value = corr 'prcorr ActiveCell.Offset(1, -1).Activate ActiveCell.Value = "STEPIT iterations" ActiveCell.Offset(0, 1).Activate If (nconv > 0) Then ActiveCell.Value = nconv Else ActiveCell.Value = ICALLS End If ActiveCell.Offset(1, -1).Activate If icode = 1 Then ActiveCell.Value = "converged" Else ActiveCell.Value = "** not converged **" End If ndf = (mc1 * mc2) - 1 - NEST 'C 'C adjust degrees of freedom for meta-cells 'C XLOGL = 0# CHISQP = 0# SUME = 0# nmeta2 = 0 For i = 1 To NCELLS SUME = SUME + FE(i) If (mask2(i) <> 0) Then nmeta2 = nmeta2 + 1 Else CHISQP = CHISQP + ((FO(i) - FE(i)) ^ 2) / FE(i) XLOGL = XLOGL + FO(i) * Log(FE(i) / total) End If Next 'C 'C now include meta-cell frequencies 'C If (nmeta > 0) Then For i = 1 To nmeta 'C next line ok? If (GE(i) <= 0#) Then GE(i) = 1E-30 CHISQP = CHISQP + ((go1(i) - GE(i)) ^ 2) / GE(i) XLOGL = XLOGL + go1(i) * Log(GE(i) / total) Next End If ndf = ndf - (nmeta2 - nmeta) If ndf > 0 Then pvalue = Application.ChiDist(CHISQL, ndf) ActiveCell.Offset(1, 0).Activate ActiveCell.Value = "Model fit " ActiveCell.Offset(1, 0).Activate ActiveCell.Value = "Likelihood ratio chi-square " ActiveCell.Offset(0, 1).Activate ActiveCell.Value = CHISQL ActiveCell.Offset(1, -1).Activate ActiveCell.Value = "p-value " ActiveCell.Offset(0, 1).Activate ActiveCell.Value = pvalue ActiveCell.Offset(1, -1).Activate ActiveCell.Value = "Peason chi-square " ActiveCell.Offset(0, 1).Activate ActiveCell.Value = CHISQP ActiveCell.Offset(1, -1).Activate ActiveCell.Value = "p-value " ActiveCell.Offset(0, 1).Activate ActiveCell.Value = Application.ChiDist(CHISQP, ndf) ActiveCell.Offset(1, -1).Activate ActiveCell.Value = "Degrees of freedom" ActiveCell.Offset(0, 1).Activate ActiveCell.Value = ndf ActiveCell.Offset(1, -1).Activate ActiveCell.Value = "Number of estimated parameters " ActiveCell.Offset(0, 1).Activate ActiveCell.Value = NEST ActiveCell.Offset(0, -1).Activate End If ActiveCell.Offset(1, 0).Activate ActiveCell.Value = "Options " ActiveCell.Offset(1, 0).Activate If (IOP(5) <> 1) Then ActiveCell.Value = "Maximum likelihood estimation " Else ActiveCell.Value = "Minimum X-squared estimation " End If ActiveCell.Offset(1, 0).Activate If (IOP(7) <> 1) Then ActiveCell.Value = "Joint estimation " Else ActiveCell.Value = "Two-step estimation " End If ActiveCell.Offset(1, 0).Activate Select Case (IOP(4)) Case 1: ActiveCell.Value = "execel(Normsdist)" '"High precision cdf (NORMP) " Case 2: ActiveCell.Value = "High precision cdf (NPROB) " Case Else: ActiveCell.Value = "Regular precision cdf(ALNORM) " End Select ActiveCell.Offset(1, 0).Activate ActiveCell.Value = mc1 ActiveCell.Offset(0, 1).Activate ActiveCell.Value = mc2 End Sub Private Sub main() 'Call init0 'C 'C Input observed data /////////////////////////////// 'C Call DATAIN 'C 'C Initialize arrays and define start values 'C Call INIT(NFIXED) 'C 'C Initialize things for STEPIT 'C Call INIT2 'C 'C Call optimization routine 'C Call STEPIT 'C 'C check flag to see if std errors could not be calculated 'C If (iflagse = 1) Then IOP(12) = 1 'C 'C Convergence reached 'C ' 'C 'C Write OUTPUT 'C 'C CALL RESCALE ' Call REPORT(1) 'C 'C Write starting values 'C ' Call STARTV(1) 'MsgBox (A(1)) End Sub Private Sub INIT2() 'C----------------------------------------------------------------------- 'C Initialize things for STEPIT 'C 'C Revised: Change(s) 'C -------- ---------------------------------------------------------- 'C 07/06/00 First written 'C----------------------------------------------------------------------- 'C 'C See STEPIT documentation 'C ICALLS = 0 NTRACE = 0 MATRIX = 105 If (IOP(12) = 1) Then MATRIX = 151 ' If (IOP(7) = 1) Then MATRIX = 121 '////// 'C write (*,*) iop(12), matrix 'C if (i = i) stop For i = 1 To NV XMAX(i) = 0# XMIN(i) = 0# DELMIN(i) = 0# DELTAX(i) = x(i) * (0.01) MASK(i) = 0 Next i 'C If nonparametric distribution, constrain latent class probabilities; 'C only xmin is important here; apparently for xmin to work, xmax 'C must have a value--so choose a big one If (IOP(11) = 0) Then For i = 1 To LC XMIN(i) = 0# XMAX(i) = 100000000# Next End If 'C 'C Constrain rho parameter to 0-1 range 'C XMIN(NEP + 1) = 0.000001 XMAX(NEP + 1) = 0.999999 'C XMAX(nep+1) = .99 End Sub Private Sub STEPIT() 'C----------------------------------------------------------------------- 'C 'C SUBROUTINE STEPIT -- GENERALIZED MAXIMIZATION/MINIMIZATION ROUTINE 'C 'C----------------------------------------------------------------------- 'C Original header comments: 'C 'C Author: J. P. Chandler 'C Physics Department, Indiana University 'C Version: 5.2 'C Written: August, 1967 'C 'C Available from Quantum Chemistry Program Exchange, I. U. 'C Chemistry Department, Bloomington, Indiana 'C 'C (Modififed for FORTRAN G OS/360 MVT) 'C----------------------------------------------------------------------- 'C Comments added by John Uebersax: 'C 'C John P. Chandler has given permission for noncommercial use of 'C the STEPIT subroutine, but requests appropriate citation. 'C The relevant reference is: 'C 'C Chandler, J. P. (1969). STEPIT--Finds local minima of a smooth 'C function of several parameters. Computer program abstract. 'C Behavioral Science, 14, 81-82. 'C 'C A user manual is available from the source mentioned in the 'C original header. 'C 'C I have tried to make only minimal changes to the original code. 'C Anything in lowercase letters is a change I have made. 'C----------------------------------------------------------------------- 'C Modified for Lahey f77L FORTRAN -- John Uebersax -- 2/10/94 'C Modified for Lahey f77P FORTRAN -- John Uebersax -- 3/14/96 'C (20 variables) 'C----------------------------------------------------------------------- Dim A, P, nt, nflag Dim VEC(70), TRIAL(70), XSAVE(70), chi(70), dx(70), SECOND(2, 2) Dim OLDVEC(70), SALVO(70), XOSC(70, 15), CHIOSC(15), NFLAT(70) 'C 'C NVMAX=50 ! f77l mod 'C KW=61 ! f77l mod MOSQUE = 15 RATIO = 10# COLIN = 0.99 NCOMP = 5 ACK = 2# 'C SIGNIF = 2.E8 ! f77l mod 'C HUGE = 1.E65 ! f77l mod SIGNIF = 2000000# ' ! f77l mod HUGE = 1E+65 ' ! f77l mod 'C JVARY = 0 40 If (NV <= 0) Then GoTo 290 ',290,50 50 NACTIV = 0 For i = 1 To NV If (MASK(i) <> 0) Then GoTo 150 ',60,150 60 If (SIGNIF * Abs(DELTAX(i)) - Abs(x(i))) > 0 Then GoTo 100 ' 70,70,100 70 If (x(i) <> 0) Then GoTo 90 ',80,90 80 DELTAX(i) = 0.01 GoTo 100 90 DELTAX(i) = 0.01 * x(i) 100 If (DELMIN(i) <> 0) Then GoTo 120 ',110,120 110 DELMIN(i) = DELTAX(i) / SIGNIF 120 If (XMAX(i) - XMIN(i) > 0) Then GoTo 140 '130,130,140 130 XMAX(i) = HUGE XMIN(i) = -HUGE 140 NACTIV = NACTIV + 1 VACD = Application.Min(XMAX(i), x(i)) x(i) = Application.Max(XMIN(i), VACD) 150 Next COMPAR = 0# 'IF(NACTIV-1)160,190,180 ee = NACTIV - 1 Select Case ee Case Is < 0: GoTo 160 Case Is = 0: GoTo 190 Case Is > 0: GoTo 180 End Select 160 For j = 1 To NV MASK(j) = 0 Next j GoTo 50 180 A = NACTIV SUB1 = 2# / (A - 1#) P = 2# * (1# / Sqr(A) / (1# - 0.5 ^ SUB1) - 1#) COMPAR = Application.Min(0.999, Abs((1# - (1# - COLIN) ^ SUB1) * (1# + P * (1# - COLIN)))) 190 If (NTRACE < 0) Then GoTo 280 ',200,200 200 'WRITE (6,210) ' 210 FORMAT(1X,' ') 'C 210 FORMAT(1X,'ENTER SUBROUTINE STEPIT. COPYRIGHT 1965 J. P. CHANDLER 'C 1, INDIANA UNIVERSITY',//,1X,'INITIAL VALUES....'/) 'C WRITE(6,220)(MASK(J),J=1,NV) 'C 220 FORMAT(/(I6,6X)/(4X,10I12)) 'C WRITE(6,230)(X(J),J=1,NV) 'C 230 FORMAT(/F10.4/(10X,10F10.4)) 'C WRITE(6,240)(XMAX(J),J=1,NV) 'C 240 FORMAT(/F10.4/(10X,10F10.4)) 'C WRITE(6,250)(XMIN(J),J=1,NV) 'C 250 FORMAT(/F10.4/(10X,10F10.4)) 'C WRITE(6,260)(DELTAX(J),J=1,NV) 'C 260 FORMAT(/E12.4/(10X,10F10.4)) 'C WRITE(6,270)(DELMIN(J),J=1,NV) 'C 270 FORMAT(/E15.4/(10X,10F10.4)) 280 Call OBJFUN '///// NF = 1 JOCK = 1 290 If (NTRACE < 0) Then GoTo 320 ',300,300 'C 'C XXXMOD JULY 1989 'C 300 i = i 'C 'C 300 WRITE(6,310)NV,NACTIV,MATRIX,NCOMP,RATIO,ACK,COLIN,COMPAR,CHISQ 'C 310 FORMAT(/1X,I3,I3,I4, 'C 1I2/1X,F5.1,10X,F5.1,F6.3,10X, 'C 2F6.3/1X,F10.8) 320 If (NV <= 0) Then GoTo 2150 ',2150,330 330 If (NTRACE <= 0) Then GoTo 360 ',360,340 'C 'C XXXMOD JULY 1989 'C 340 i = i 'C 'C 340 WRITE(6,350) 'C 350 FORMAT(/10(1X,1H*)) 'C 360 For i = 1 To NV dx(i) = DELTAX(i) VEC(i) = 0# For j = 1 To 20 370 ERR(i, j) = 0# Next j Next i CHIOLD = CHISQ NOSC = 0 'C 380 NCIRC = 0 NZIP = 0 'C 'C MAIN DO LOOP FOR CYCLING THROUGH THE VARIABLES. 'C FIRST TRIAL STEP WITH EACH VARIABLE IS SEPARATE. 390 NACK = 0 For i = 1 To NV OLDVEC(i) = VEC(i) VEC(i) = 0# TRIAL(i) = 0# If (MASK(i) = 0) Then GoTo 410 ' 400,410,400 400 VEC(i) = -0# NFLAT(i) = 1 GoTo 1350 410 NACK = NACK + 1 Save = x(i) If (SIGNIF * Abs(dx(i)) - Abs(x(i))) <= 0 Then GoTo 580 ',580,420 420 x(i) = Save + dx(i) JVARY = 0 If (JOCK <= 0) Then GoTo 440 ',440,430 430 JOCK = 0 JVARY = i 440 nflag = 1 If (x(i) - XMIN(i)) <= 0 Then GoTo 460 ',460,450 450 If (x(i) - XMAX(i)) < 0 Then GoTo 470 ',460,460 460 nflag = nflag + 3 GoTo 490 470 Call OBJFUN '/// If (ICALLS <= IMAX) Then GoTo 10 Exit Sub 10 'next NF = NF + 1 JVARY = i CHIME = CHISQ ' IF(CHISQ-CHIOLD)620,480,490 ee = CHISQ - CHIOLD Select Case ee Case Is < 0: GoTo 620 Case Is = 0: GoTo 480 Case Is > 0: GoTo 490 End Select 480 nflag = nflag + 1 490 x(i) = Save - dx(i) If (x(i) - XMIN(i)) <= 0 Then GoTo 590 ',590,500 500 If (x(i) - XMAX(i)) >= 0 Then GoTo 590 '510,590,590 510 Call OBJFUN '// NF = NF + 1 JVARY = i 'IF(CHISQ-CHIOLD)610,520,530 ee = CHISQ - CHIOLD Select Case ee Case Is < 0: GoTo 610 Case Is = 0: GoTo 520 Case Is > 0: GoTo 530 End Select 520 nflag = nflag + 1 530 'IF(NFLAG-3)540,580,590 ee = nflag - 3 Select Case ee Case Is < 0: GoTo 540 Case Is = 0: GoTo 580 Case Is > 0: GoTo 590 End Select 540 If ((CHISQ - CHIME) * (CHIME - 2# * CHIOLD + CHISQ)) = 0 Then GoTo 590 '550,590,550 550 TRIAL(i) = 0.5 * dx(i) * (CHISQ - CHIME) / (CHIME - 2# * CHIOLD + CHISQ) VEC(i) = TRIAL(i) / Abs(dx(i)) NFLAT(i) = 0 x(i) = Save + TRIAL(i) Call OBJFUN '/// NF = NF + 1 If (CHISQ - CHIOLD) >= 0 Then GoTo 570 '560,570,570 560 CHIOLD = CHISQ JOCK = 1 GoTo 600 570 TRIAL(i) = 0# VEC(i) = 0# GoTo 590 580 VEC(i) = -0# NFLAT(i) = 1 590 x(i) = Save 600 NCIRC = NCIRC + 1 ' If (NCIRC - NACTIV) >= 0 Then GoTo 1430 '690,1430,1430 ee = NCIRC - NACTIV Select Case ee Case Is < 0: GoTo 690 Case Is = 0: GoTo 1430 Case Is > 0: GoTo 1430 End Select 610 dx(i) = -dx(i) 'C 'C A LOWER VALUE HAS BEEN FOUND. HENCE THIS VARIABLE WILL CHANGE. 'C 620 NCIRC = 0 DEL = dx(i) 630 CHIME = CHIOLD CHIOLD = CHISQ VEC(i) = VEC(i) + DEL / Abs(dx(i)) NFLAT(i) = 0 TRIAL(i) = TRIAL(i) + DEL DEL = ACK * DEL Save = x(i) x(i) = Save + DEL If (x(i) - XMIN(i)) <= 0 Then GoTo 680 ',680,640 640 If (x(i) - XMAX(i)) >= 0 Then GoTo 680 '650,680,680 650 Call OBJFUN '////////////// NF = NF + 1 If (CHISQ - CHIOLD) < 0 Then GoTo 630 ',660,660 660 ZZZ = ACK * CHIME - (ACK + 1#) * CHIOLD + CHISQ If (ZZZ = 0#) Then GoTo 661 CINDER = (0.5 / ACK) * (ACK ^ 2 * CHIME - (ACK ^ 2 - 1#) * CHIOLD - CHISQ) / ZZZ GoTo 662 661 CINDER = 0# 662 x(i) = Save + CINDER * DEL Call OBJFUN '/////// NF = NF + 1 If (CHISQ - CHIOLD) >= 0 Then GoTo 680 '670,680,680 670 CHIOLD = CHISQ TRIAL(i) = TRIAL(i) + CINDER * DEL VEC(i) = VEC(i) + CINDER * DEL / Abs(dx(i)) GoTo 690 680 x(i) = Save 690 If (NZIP - 1) < 0 Then GoTo 1340 ',700,700 700 If (Abs(VEC(i)) - ACK) < 0 Then GoTo 750 ',710,710 710 dx(i) = ACK * Abs(dx(i)) VEC(i) = VEC(i) / ACK OLDVEC(i) = OLDVEC(i) / ACK For j = 1 To MOSQUE 720 ERR(i, j) = ERR(i, j) / ACK Next j If (NTRACE) > 0 Then '750,750,730 730 'WRITE(6,740)I,DX(I) ' 740 FORMAT(1X,I3,F10.5) End If 750 SUMO = 0# SUMV = 0# For j = 1 To NV SUMO = SUMO + OLDVEC(j) ^ 2 760 SUMV = SUMV + VEC(j) ^ 2 Next j If (SUMO * SUMV) <= 0 Then GoTo 1340 ',1340,770 770 SUMO = Sqr(SUMO) SUMV = Sqr(SUMV) COSINE = 0# For j = 1 To NV COSINE = COSINE + (OLDVEC(j) / SUMO) * (VEC(j) / SUMV) Next j 'IF(NZIP-1)1340,790,800 ee = NZIP - 1 Select Case ee Case Is < 0: GoTo 1340 Case Is = 0: GoTo 790 Case Is > 0: GoTo 800 End Select 790 'If (NACK - NACTIV) < 0 Then GoTo 1340 ',820,820 ee = NACK - NACTIV Select Case ee Case Is < 0: GoTo 1340 Case Is = 0: GoTo 820 Case Is > 0: GoTo 820 End Select 800 If (NACK - NACTIV) < 0 Then GoTo 820 ',810,810 810 If (NZIP - NCOMP) >= 0 Then GoTo 830 ' 820,830,830 820 If (COSINE - COMPAR) < 0 Then GoTo 1340 ',830,830 'C 'C SIMON SAYS, TAKE AS MANY GIANT STEPS AS POSSIBLE... 'C 830 If (NTRACE) > 0 Then '860,860,840 840 'WRITE(6,850)CHIOLD,(VEC(J),J=1 to I) ' 850 FORMAT(1X,F10.8,5X,10F9.2/,1X,10F9.2) End If 860 NGIANT = 0 NTRY = 0 NRETRY = 0 KL = 1 NOSC = NOSC + 1 If (NOSC - MOSQUE) <= 0 Then GoTo 890 ',890,870 870 NOSC = MOSQUE For k = 2 To MOSQUE CHIOSC(k - 1) = CHIOSC(k) For j = 1 To NV XOSC(j, k - 1) = XOSC(j, k) ERR(j, k - 1) = ERR(j, k) Next j Next k 890 For j = 1 To NV XOSC(j, NOSC) = x(j) ERR(j, NOSC) = VEC(j) / SUMV Next j CHIOSC(NOSC) = CHIOLD If (NOSC - 3) < 0 Then GoTo 960 ',910,910 'C 'C SEARCH FOR A PREVIOUS SUCCESSFUL GIANT STEP IN A DIRECTION MORE 'C NEARLY PARALLEL TO THE DIRECTION OF THE PROPOSED STEP THAN WAS THE 'C IMMEDIATELY PREVIOUS ONE. 'C 910 COXCOM = 0# For j = 1 To NV COXCOM = COXCOM + ERR(j, NOSC) * ERR(j, NOSC - 1) Next j NAH = NOSC - 2 930 NTRY = 0 For k = KL To NAH NRETRY = NAH - k COSINE = 0# For j = 1 To NV COSINE = COSINE + ERR(j, NOSC) * ERR(j, k) Next j If (COSINE - COXCOM) > 0 Then GoTo 970 '950,950,970 950 Next k 960 CHIBAK = chi(i) GoTo 1020 970 NTRY = 1 KL = k + 1 If (NTRACE) <= 0 Then GoTo 1000 ',1000,980 980 nt = NOSC - k ' WRITE(6,990)NT ' 990 FORMAT(/1X,I2) 1000 For j = 1 To NV SALVO(j) = TRIAL(j) TRIAL(j) = (x(j) - XOSC(j, k)) / ACK Next j CHIBAK = CHIOLD + (CHIOSC(k) - CHIOLD) / ACK 'C 1020 For j = 1 To NV XSAVE(j) = x(j) TRIAL(j) = ACK * TRIAL(j) If (MASK(j)) = 0 Then '1040,1030,1040 TRAAL = x(j) + TRIAL(j) TRAIL = Application.Min(TRAAL, XMAX(j)) x(j) = Application.Max(TRAIL, XMIN(j)) End If Next JOCK = 0 JVARY = 0 Call OBJFUN '//////// NF = NF + 1 If (CHISQ - CHIOLD) >= 0 Then GoTo 1080 '1050,1080,1080 1050 CHIBAK = CHIOLD CHIOLD = CHISQ NGIANT = NGIANT + 1 If (NTRACE) <= 0 Then GoTo 1020 ',1020,1060 1060 ' WRITE(6,1070)CHISQ,(X(J),J=1,NV) ' 1070 FORMAT(1X,F10.8,/(10(1X,F10.5))) GoTo 1020 'C 1080 If (NRETRY) <= 0 Then GoTo 1100 ',1100,1090 1090 If (NGIANT) <= 0 Then GoTo 1150 ',1150,1100 1100 CINDER = (0.5 / ACK) * (ACK ^ 2 * CHIBAK - (ACK ^ 2 - 1#) * CHIOLD - CHISQ) / (ACK * CHIBAK - (ACK + 1#) * CHIOLD + CHISQ) For j = 1 To NV If (MASK(j) = 0) Then '1120,1110,1120 CANDER = XSAVE(j) + CINDER * TRIAL(j) XIJ = Application.Min(CANDER, XMAX(j)) x(j) = Application.Max(XIJ, XMIN(j)) End If Next JOCK = 0 JVARY = 0 Call OBJFUN '///// NF = NF + 1 If (CHISQ - CHIOLD) < 0 Then GoTo 1280 ',1130,1130 1130 If (NGIANT <> 0) Then GoTo 1170 ',1140,1170 1140 If (NTRY = 0) Then GoTo 1170 '1150,1170,1150 1150 For j = 1 To NV TRIAL(j) = SALVO(j) x(j) = XSAVE(j) Next j GoTo 1190 1170 For j = 1 To NV TRIAL(j) = TRIAL(j) / ACK x(j) = XSAVE(j) Next j 1190 If (NTRACE <= 0) Then GoTo 1240 ',1240,1200 ' 1200 WRITE(6,1210)CHIOLD,NGIANT ' 1210 FORMAT(1X,F10.8) ' WRITE(6,1220)(X(J),J=1 to NV) ' 1220 FORMAT((10(1X,F10.5))) ' WRITE(6,1230) ' 1230 FORMAT(/1X,' ') 1240 If (NGIANT > 0) Then GoTo 1310 '1250,1250,1310 1250 If (NRETRY > 0) Then GoTo 930 '1260,1260,930 1260 If (NTRY = 0) Then GoTo 1330 '1270,1330,1270 1270 NTRY = 0 GoTo 960 'C 1280 CHIOLD = CHISQ JOCK = 1 If (NTRACE <= 0) Then GoTo 1310 ',1310,1290 1290 STEPS = (NGIANT) + CINDER ' WRITE(6,1300)CHIOLD,STEPS '1300 FORMAT(1X,F10.8) ' WRITE(6,1220)(X(J),J=1 to NV) ' WRITE(6,1230) 1310 If (NTRY = 0) Then GoTo 380 ' 1320,380,1320 1320 NOSC = 0 GoTo 380 1330 NOSC = Application.Max(NOSC - 1, 0) 1340 chi(i) = CHIOLD 1350 Next 'C 'C ANOTHER CYCLE THROUGH THE VARIABLES HAS BEEN COMPLETED. 'C PRINT ANOTHER LINE OF TRACES. 'C If (NTRACE > 0) Then '1370,1370,1360 1360 'WRITE(6,850)CHIOLD,(VEC(J),J=1 to NV) End If 1370 ' Next 1380 If (NZIP <> 0) Then GoTo 1420 ',1390,1420 1390 If (NTRACE <= 0) Then GoTo 1420 ',1420,1400 1400 'WRITE(6,1220)(X(J),J=1 to NV) 'WRITE(6,1410) '1410 FORMAT(1H ) 1420 NZIP = NZIP + 1 GoTo 390 'C 'C A MINIMUM HAS BEEN FOUND. PRINT THE REMAINING TRACES. 'C 1430 If (NTRACE <= 0) Then GoTo 1450 ' ,1450,1440 1440 'WRITE(6,850)CHIOLD,(VEC(J),J=1 to I) 1450 If (NTRACE <= 0) Then GoTo 1470 ',1470,1460 1460 'WRITE(6,1220)(X(J),J=1 to NV) 'WRITE(6,1230) 'C 'C DECREASE THE SIZE OF THE STEPS FOR ALL VARIABLES. 'C 1470 'next 1480 NOSC = 0 NGATE = 1 For j = 1 To NV If (MASK(j) <> 0) Then GoTo 1520 ',1490,1520 1490 If (NFLAT(j) > 0) Then GoTo 1520 '1500,1500,1520 1500 If (Abs(dx(j)) - Abs(DELMIN(j))) <= 0 Then GoTo 1520 ',1520,1510 1510 NGATE = 0 1520 dx(j) = dx(j) / RATIO Next If (NGATE > 0) Then GoTo 1600 '1530,1530,1600 1530 If (NTRACE <= 0) Then GoTo 1570 ',1570,1540 1540 'WRITE(6,1550)(DX(J),J=1,NV) '1550 FORMAT(10(1X,1H*)/,(10(1X,F10.5))) ' WRITE(6,1560) '1560 FORMAT(//) 1570 GoTo 380 'C1580 WRITE(6,1590)(DX(J),J=1,NV) ! f77l mod 'C1590 FORMAT (1X,'OPERAT. TERM.'//10(1X,F10.5)) ! f77l mod 1600 CHISQ = CHIOLD If (NTRACE < 0) Then GoTo 1630 ',1610,1610 1610 'WRITE(6,1620)NF nconv = NF icode = 1 '/// ' write (*,*) 'Convergence reached' ICALLS = 0 ' 1620 FORMAT(/1X,I5,23H FUNCTION COMPUTATIONS ) 1630 ' next 1640 If (Int(Abs(MATRIX - 100)) - 50) > 0 Then GoTo 2160 ' 1650,1650,2160 1650 If (NACTIV - NV) <> 0 Then GoTo 2160 ',1660,2160 'C 'C COMPUTE THE STANDARD ERRORS AND THE CORRELATIONS. 'C 1660 FAC = RATIO ^ (MATRIX - 100) ' write (*, *) 'Standard errors: information matrix' For i = 1 To NV dx(i) = Abs(FAC * dx(i)) XSAVE(i) = x(i) JVARY = 0 For j = 1 To 2 x(i) = XSAVE(i) + dx(i) Call OBJFUN '//// NF = NF + 1 JVARY = i SECOND(1, j) = CHISQ dx(i) = -dx(i) Next j x(i) = XSAVE(i) ERR(i, i) = (SECOND(1, 1) - 2# * CHIOLD + SECOND(1, 2)) / dx(i) ^ 2 Next i For i = 2 To NV im = i - 1 For j = 1 To im For k = 1 To 2 x(i) = XSAVE(i) + dx(i) JVARY = 0 For l = 1 To 2 x(j) = XSAVE(j) + dx(j) Call OBJFUN '///// NF = NF + 1 JVARY = j SECOND(k, l) = CHISQ x(j) = XSAVE(j) dx(j) = -dx(j) Next l x(i) = XSAVE(i) dx(i) = -dx(i) Next k ERR(i, j) = 0.25 * (SECOND(1, 1) - SECOND(1, 2) - SECOND(2, 1) + SECOND(2, 2)) / Abs(dx(i) * dx(j)) ERR(j, i) = ERR(i, j) Next j Next i If (NTRACE < 0) Then GoTo 1770 ',1720,1720 1720 'WRITE(6,1730) '1730 FORMAT(40HSIZES OF INCREMENTS TO BE USED BELOW....) ' WRITE(6,1740)(DX(J),J=1,NV) '1740 FORMAT(/(10(1X,F10.5))) ' WRITE(6,1750) '1750 FORMAT(/45H MATRIX OF THE SECOND PARTIAL DERIVATIVES.... /) For i = 1 To NV 1760 ' WRITE(6,1740)(ERR(I,J),J=1,I) Next 1770 For i = 1 To NV For j = 1 To i '/// If (ERR(i, j) = 0) Then GoTo 1790 '1780,1790,1780 Next j, i GoTo 1810 1790 ' WRITE(6,1800) '1800 FORMAT(////123H MsgBox ("THE ABOVE MATRIX CONTAINS ONE OR MORE ZEROES. A LARGER VALUE OF 'MATRIX' SHOULD BE TRIED, TO SEE IF THEY ARE LEGITIMATE.") iflagse = 1 '//// Exit Sub GoTo 2160 '//// 'C 'C INVERT THE MATRIX USING SYMINV2 (COMM. OF THE A.C.M. 6, P. 67). 'C 1810 det = 1# ' write (*, *) 'Standard errors: inverting matrix' DETLOG = 0# For j = 1 To NV SALVO(j) = 1# Next j For i = 1 To NV 'for 1970 I=1 to NV ??????? BIGAJJ = 0# For j = 1 To NV If (SALVO(j) = 0) Then GoTo 1850 '1830,1850,1830 1830 If (Abs(ERR(j, j)) - BIGAJJ) <= 0 Then GoTo 1850 ',1850,1840 1840 BIGAJJ = Abs(ERR(j, j)) k = j 1850 Next If (BIGAJJ <> 0) Then GoTo 1870 ',1860,1870 1860 det = 0# GoTo 1980 1870 SALVO(k) = 0# 'C IF (DET > 1.E48 or DET < 1.E-48) GO TO 1871 ! f77l mod If (det > 1E+48 Or det < 1E-48) Then 'GO TO 1871 ! f77l mod ' GO TO 1875 1871 ' WRITE (6,1872) 'C1872 FORMAT (///,'DETERMINANT COMPUTATION WILL CAUSE OVERFLOW.') ' write (2, 1872) '1872 FORMAT(//, MsgBox ("Problem: Calculation of the determinant of the matrix of second-derivatives would cause overflow. This could mean that the model is weakly identified. Standard error estimation is suppressed.") 'C 'C flag: standard error calculated is suppressed 'C iflagse = 1 GoTo 2160 End If 1875 det = det * ERR(k, k) DETLOG = DETLOG + Log(Abs(ERR(k, k))) / 2.303 TRIAL(k) = 1# / ERR(k, k) ERR(k, k) = 0# XSAVE(k) = 1# M = k - 1 If (M <= 0) Then GoTo 1910 ',1910,1880 1880 For j = 1 To M XSAVE(j) = ERR(k, j) TRIAL(j) = ERR(k, j) * TRIAL(k) ' IF(SALVO(J))1860,1900,1890 Select Case SALVO(j) Case Is < 0: GoTo 1860 Case Is = 0: GoTo 1900 Case Is > 0: GoTo 1890 End Select 1890 TRIAL(j) = -TRIAL(j) 1900 ERR(k, j) = 0# Next 1910 M = k + 1 If (M - NV) > 0 Then GoTo 1960 '1920,1920,1960 1920 For j = M To NV XSAVE(j) = ERR(j, k) 'IF(SALVO(J))1860,1930,1940 Select Case SALVO(j) Case Is < 0: GoTo 1860 Case Is = 0: GoTo 1930 Case Is > 0: GoTo 1940 End Select 1930 XSAVE(j) = -XSAVE(j) 1940 TRIAL(j) = -ERR(j, k) * TRIAL(k) ERR(j, k) = 0# Next j 1960 For j = 1 To NV For k = j To NV 1970 ERR(k, j) = ERR(k, j) + XSAVE(j) * TRIAL(k) '??????? Next k Next j Next i 'IF(DET)2000,1980,2020 Select Case det Case Is < 0: GoTo 2000 Case Is = 0: GoTo 1980 Case Is > 0: GoTo 2020 End Select 1980 'WRITE(6,1990) 1990 MsgBox ("ERROR MATRIX IS SINGULAR. 'MATRIX' SHOULD PROBABLY BE INCREASED.") GoTo 2150 2000 'WRITE(6,2010) 2010 MsgBox ("ERROR MATRIX IS NEGATIVE DEFINITE. MATRIX SHOULD PROBABLY BE DECREASED.") 2020 If (NTRACE < 0) Then GoTo 2050 ',2030,2030 2030 'WRITE(6,2040)DET,DETLOG '2040 FORMAT(/1X,F10.5,10X,F10.5) 2050 For i = 1 To NV For j = 1 To i ERR(i, j) = ERR(i, j) * 2# ERR(j, i) = ERR(i, j) Next j If (ERR(i, i) > 0) Then GoTo 2090 '2070 ',2070,2090 2070 'WRITE(6,2080)ERR(I,I) '2080 FORMAT(///50H NEGATIVE OR ZERO MEAN SQUARE ERROR ENCOUNTERED..., ' 13X,F10.8/39H 'MATRIX' SHOULD PROBABLY BE DECREASED. ///) 2090 XSAVE(i) = Sgn(ERR(i, i)) * Sqr(Abs(ERR(i, i))) 'DSIGN(sqr(abs(ERR(I,I))),ERR(I,I)) ??? Next i If (NTRACE < 0) Then GoTo 2160 ',2100,2100 2100 'WRITE(6,2110) '2110 FORMAT(/////20H STANDARD ERRORS.... ) 'WRITE(6,1740)(XSAVE(J),J=1,NV) 'C 'C write standard errors to output file 'C 'C WRITE(2,2110) 'C WRITE(2,1740)(XSAVE(J),J=1,NV) ' WRITE(6,2120) 2120 ' FORMAT(/////45H LOWER TRIANGLE OF THE CORRELATION MATRIX.... /) ' for 2140 I=2,NV ' IM=I-1 ' for 2130 J=1 to IM '2130 TRIAL(J)=ERR(I,J)/abs(XSAVE(I)*XSAVE(J)) '2140 WRITE(6,1740)(TRIAL(J),J=1 to IM) 2150 'WRITE(6,1620)NF 'C 2160 ' 2190 JVARY = 0 Call OBJFUN '///// If (NTRACE >= 0) Then ' goto 2230 ',2200,2200 2200 'WRITE(6,2210)(X(J),J=1,NV) '2210 FORMAT(/1X,/(7(1X,E16.9))) ' 'WRITE(6,2220)CHISQ '2220 'FORMAT(/1X,F10.8) End If ' 2230 End Sub Private Sub OBJFUN() 'C----------------------------------------------------------------------- 'C SUBROUTINE OBJFUN -- Objective function to minimize/maximize, in this 'C case, minimize likelihood ratio chi-square, which is the same as 'C maximizing the log-likelihood 'C 'C Revised: Change(s) 'C -------- ---------------------------------------------------------- 'C 04/24/92 Improvements to code 'C 07/15/00 Warn on negative G-squared 'C----------------------------------------------------------------------- 'C Unload parameter vector 'C INDX = 0 'C 'C Latent class prevalences, if nonparametric distribution 'C If (IOP(11) = 0) Then RTOT = 0# For i = 1 To LC INDX = INDX + 1 r(i) = x(INDX) RTOT = RTOT + r(i) Next i 'C 'C rescale prevalences so that they sum to 1.0 'C For i = 1 To LC r(i) = r(i) / RTOT Next i End If 'C 'C Rho 'C For i = 1 To NDIS INDX = INDX + 1 A(i) = x(INDX) Next i 'C 'C Thresholds, if joint estimation 'C 18 If (IOP(7) <> 1) Then For j = 2 To mc1 INDX = INDX + 1 t(1, j) = x(INDX) 'C ' write (*,*) t(1,j) Next j For j = 2 To mc2 INDX = INDX + 1 t(2, j) = x(INDX) 'C ' write (*,*) t(2,j) 21 Next j End If 'C 'C Calculate conditional probabilities 'C 'C 'C Variable 1 'C i = 1 For j = 1 To LC 'C 'C no bias model in polycorr 'C OFS = 0# 'C 'C sd is measurement error standard deviation 'C pcorr = A(1) aa2 = (1# - pcorr) / pcorr SD = aa2 ^ 0.5 'C 'C if IT model, use same thresholds for both raters 'C If (IOP(6) = 1) Then ii = 1 Else ii = i End If 'C 'C manifest category 1 'C 'C Z = (T(II, 2) + OFS - B(J)) / SD 'C 'C use latent response variable thresholds 'C If (nflag = 0 Or i = 1) Then TTT = t(ii, 2) / pcorr ^ 0.5 z = (TTT + OFS - B(j)) / SD P(i, j, 1) = phi(z) Else TTT = -t(ii, 2) / pcorr ^ 0.5 z = (TTT + OFS - B(j)) / SD P(i, j, 1) = 1# - phi(z) End If Sum = P(i, j, 1) 'C 'C manifest categories 2 to (mc1 - 1) 'C For k = 2 To mc1 - 1 'C Z1 = (T(II, K ) + OFS - B(J)) / SD 'C Z2 = (T(II, K+1) + OFS - B(J)) / SD 'C 'C use latent response variable thresholds 'C If (nflag = 0 Or i = 1) Then TTT1 = t(ii, k) / pcorr ^ 0.5 TTT2 = t(ii, k + 1) / pcorr ^ 0.5 Z1 = (TTT1 + OFS - B(j)) / SD Z2 = (TTT2 + OFS - B(j)) / SD P(i, j, k) = phi(Z2) - phi(Z1) Else TTT1 = -t(ii, k) / pcorr ^ 0.5 TTT2 = -t(ii, k + 1) / pcorr ^ 0.5 Z1 = (TTT1 + OFS - B(j)) / SD Z2 = (TTT2 + OFS - B(j)) / SD P(i, j, k) = phi(Z1) - phi(Z2) End If Sum = Sum + P(i, j, k) Next k 'C 'C manifest category mc1 'C P(i, j, mc1) = 1 - Sum Next j 'C 'C Variable 2 'C i = 2 For j = 1 To LC 'C 'C no bias model in polycorr 'C OFS = 0# 'C 'C sd is measurement error standard deviation 'C pcorr = A(1) aa2 = (1# - pcorr) / pcorr SD = aa2 ^ 0.5 'C 'C if IT model, use same thresholds for both raters 'C If (IOP(6) = 1) Then ii = 1 Else ii = i End If 'C 'C manifest category 1 'C 'C Z = (T(II, 2) + OFS - B(J)) / SD 'C 'C use latent response variable thresholds 'C If (nflag = 0 Or i = 1) Then TTT = t(ii, 2) / pcorr ^ 0.5 z = (TTT + OFS - B(j)) / SD P(i, j, 1) = phi(z) Else TTT = -t(ii, 2) / pcorr ^ 0.5 z = (TTT + OFS - B(j)) / SD P(i, j, 1) = 1# - phi(z) End If Sum = P(i, j, 1) 'C 'C manifest categories 2 to (mc2 - 1) 'C For k = 2 To mc2 - 1 'C Z1 = (T(II, K ) + OFS - B(J)) / SD 'C Z2 = (T(II, K+1) + OFS - B(J)) / SD 'C 'C use latent response variable thresholds 'C If (nflag = 0 Or i = 1) Then TTT1 = t(ii, k) / pcorr ^ 0.5 TTT2 = t(ii, k + 1) / pcorr ^ 0.5 Z1 = (TTT1 + OFS - B(j)) / SD Z2 = (TTT2 + OFS - B(j)) / SD P(i, j, k) = phi(Z2) - phi(Z1) Else TTT1 = -t(ii, k) / pcorr ^ 0.5 TTT2 = -t(ii, k + 1) / pcorr ^ 0.5 Z1 = (TTT1 + OFS - B(j)) / SD Z2 = (TTT2 + OFS - B(j)) / SD P(i, j, k) = phi(Z1) - phi(Z2) End If Sum = Sum + P(i, j, k) Next 'C 'C manifest category mc2 'C P(i, j, mc2) = 1 - Sum Next 'C 'C Expected frequencies 'C 'C 'C reset all meta-cell expected frequencies to 0 'C If (nmeta > 0) Then For i = 1 To nmeta GE(i) = 0# Next End If ETOTAL = 0# For i = 1 To NCELLS 'C IF (IOP(2) = 1 and FO(I) = 0.0) then go to 40 Sum = 0# For j = 1 To LC PROD = r(j) For k = 1 To NR PROD = PROD * P(k, j, IVEC(i, k)) Next k Sum = Sum + PROD Next j FE(i) = Sum * total ETOTAL = ETOTAL + FE(i) 'C 'C update meta-cell expected frequency, if appropriate 'C If (mask2(i) <> 0) Then mm = mask2(i) GE(mm) = GE(mm) + FE(i) End If Next i 'C 'C Evaluate fit 'C CHISQL = 0# CHISQP = 0# 'C XLOGL = 0. For ii = 1 To NCELLS 'C 'C skip to end of loop if cells belongs to a meta-cell 'C If (mask2(ii) <> 0) Then GoTo 50 'C some question about following?? try commenting out for now 'C is there a limit on how small the argument of dlog() can be 'C check reference or else experiment 'C IF (FE(II) < 1.D-10) FE(II) = 1.D-10 'C IF (FE(II) <= 0.) then 'C write (2, *) 'ii = ', ii 'C write (2, *) 'x(1) ', x(1) 'C call pause 'C end if If (FE(ii) <= 0#) Then FE(ii) = 1E-160 'C if (fe(ii) < 1.d-10 and fo(ii) <> 0) then 'C write (*, *) fe(ii), dlog(fe(ii)) 'C end if 'C 'C G-squared 'C If (FO(ii) <> 0#) Then CHISQL = CHISQL - 2# * FO(ii) * Log(FE(ii) / FO(ii)) 'C XLOGL = XLOGL + FO(II) * LOG(FE(II) / TOTAL) End If 'C 'C X-squared, if applicable 'C If (IOP(5) = 1) Then CHISQP = CHISQP + ((FO(ii) - FE(ii)) ^ 2) / FE(ii) End If 50 Next ii 'C 'C now consider the meta-cells 'C For ii = 1 To nmeta If (GE(ii) <= 0#) Then GE(ii) = 1E-160 If (IOP(5) = 0) Then If (go1(ii) <> 0#) Then CHISQL = CHISQL - 2# * go1(ii) * Log(GE(ii) / go1(ii)) 'C XLOGL = XLOGL + go1(II) * LOG(ge(II) / TOTAL) End If Else CHISQP = CHISQP + ((go1(ii) - GE(ii)) ^ 2) / GE(ii) End If Next ii If (CHISQL < 0#) Then ' write (*, 6200) MsgBox ("ERROR: Run terminated due to negative G-squared value. This could be due to a very small marginal proportion. Consider collapsing the table.") If (IOP(7) <> 1) Then 'write (*, 6202) ///// MsgBox (" 2-step estimation might avoid the problem.") End If ' write (2, 2200) '2200 format('ERROR: Negative G-squared value calculated',//, ' 1 'Values of parameter vector at time of error:',//, ' 2 'Parameter no. Value',/, ' 3 '------------- ------------') ' for 55 i = 1 to nv ' write (2, 2202) i, x(i) ' 55 next ' 2202 format(5x, i7, 5x, f10.5) GoTo 180 End If If (IOP(5) = 0) Then CHISQ = CHISQL Else CHISQ = CHISQP End If 'C CHISQ = -XLOGL ICALLS = ICALLS + 1 ' WRITE (*,6000) ICALLS, CHISQ If (ICALLS < IMAX) Then GoTo 200 'C 'C Maximum number of iterations reached without convergence 'C 180 MsgBox ("**Convergence not reached** " + Str(ICALLS)) 'WRITE (*, 6010) 'C 'C Write start values 'C ' CALL STARTV (1) 'C 'C Write OUTPUT 'C ' CALL REPORT (0) End 200 ' Next ' 6000 FORMAT(1X, I5, 5X, F16.8) ' 6010 FORMAT(/,' **Convergence not reached**') End Sub Private Sub init0() TITLE = "" '繰り返し最大数 IMAX = 5000 ' READ (5, 5001) mc1 ' READ (5, 5001) mc2 NR = 2 ' No. of items or raters For i = 1 To NR BIAS(i) = 0# Next i '外部ファイルからの初期値の読み込み 'C IOP(1) Use file start values (0=no, 1=yes) IOP(1) = 0 'C x IOP(2) Indexed frequencies (0=no, 1=yes) IOP(2) = 0 'C x IOP(3) Equal measurement error (0=no, 1=yes) IOP(3) = 1 If (IOP(3) = 1) Then NDIS = 1 'C IOP(4) cdf precision: 0=regular, 1=high, 2 =NPROB IOP(4) = 2 'C iop(5) Maximum-likelihood (0) or minumum X-squared (1) criterion IOP(5) = 0 'C x IOP(6) Constrained thresholds (0=no; 1 = yes) IOP(6) = 0 nt = NR If (IOP(6) > 0) Then nt = 1 'C IOP(7) fiml(joint ML estimation) (0) or 2-step (1) estimation IOP(7) = 1 'C x IOP(8) Print expected frequencies (0=no, 1=yes) IOP(8) = 1 '固定しているが,0 のほうがいいだろう 'C x IOP(9) Print conditional rating probabilities (0=no, 1=yes) IOP(9) = 0 'C x IOP(10) Print classification information (0=no, 1=yes) IOP(10) = 0 'C x IOP(11) Distribution type (0=nonpar, 1=normal) IOP(11) = 1 'C IOP(12) Suppresss std err calculation (0=no, 1=yes) IOP(12) = 0 If (IOP(5) = 1) Then IOP(12) = 1 'C IOP(13) Not used yet. IOP(13) = 0 'C IOP(14) High precision format of expected frequencies IOP(14) = 1 'C IOP(15) Not used yet. IOP(15) = 0 'C Latent trait range tlimit = 5 If (tlimit <= 0#) Then tlimit = 5# If (tlimit > 10#) Then tlimit = 10# 'C No. of latent classes/quadrature points 'C LC = 51 If (LC <= 0) Then LC = 51 If (LC > 81) Then LC = 81 'C Number of meta-calls 'C nmeta = 0 If (nmeta > 20) Then nmeta = 20 'C no fixed parameters in polycorr 'C NFIXED = 0 'C nep = number of estimated latent class prevalences 'C iop(11)=1 固定だから nep=0 If (IOP(11) = 0) Then NEP = LC Else NEP = 0 End If End Sub Private Sub INIT(NFIXED) 'C----------------------------------------------------------------------- 'C Initialize arrays, define start values 'C クロス表の作成と最初の計算 'C Revised: Change(s) 'C -------- ---------------------------------------------------------- 'C 04/23/92 Made separate subroutine 'C----------------------------------------------------------------------- NFIXED = 0 'C 'C Generate index vectors for 2-way table 'C INDX = 0 For i = 1 To mc1 For j = 1 To mc2 INDX = INDX + 1 IVEC(INDX, 1) = i IVEC(INDX, 2) = j Next j Next i 'C 'C Calculate 1-way marginals and cumulative marginals 'C For i = 1 To mc1 FMARG(1, i) = 0# CFMARG(1, i) = 0# Next i For i = 1 To mc2 FMARG(2, i) = 0# CFMARG(2, i) = 0# Next i For l = 1 To NCELLS ii = IVEC(l, 1) jj = IVEC(l, 2) FMARG(1, ii) = FMARG(1, ii) + FO(l) FMARG(2, jj) = FMARG(2, jj) + FO(l) Next l CFMARG(1, 1) = FMARG(1, 1) CFMARG(2, 1) = FMARG(2, 1) For i = 2 To mc1 CFMARG(1, i) = CFMARG(1, i - 1) + FMARG(1, i) Next For i = 2 To mc2 CFMARG(2, i) = CFMARG(2, i - 1) + FMARG(2, i) Next 'C 'C make sure both items/raters use all levels 'C For i = 1 To mc1 If (FMARG(1, i) = 0#) Then MsgBox (" ERROR: All rating levels were not used by both items/raters. Polychoric correlation cannot be calculated. See documentation for details.") End Exit Sub End If Next For i = 1 To mc2 If (FMARG(2, i) = 0#) Then MsgBox (" ERROR: All rating levels were not used by both items/raters. Polychoric correlation cannot be calculated. See documentation for details.") End Exit Sub End If Next i 'C 'C Calculate G-squared and X-squared for independence model 'C For i = 1 To mc1 For j = 1 To mc2 For k = 1 To NCELLS If ((IVEC(k, 1) = i) And (IVEC(k, 2) = j)) Then FE(k) = FMARG(1, i) * FMARG(2, j) / total End If Next k Next j Next i chisql0 = 0# chisqp0 = 0# For i = 1 To NCELLS If (FE(i) < 0.0000000001) Then FE(i) = 0.0000000001 chisqp0 = chisqp0 + ((FO(i) - FE(i)) ^ 2) / FE(i) If (FO(i) <> 0#) Then chisql0 = chisql0 - 2# * FO(i) * Log(FE(i) / FO(i)) End If Next i ndf0 = (mc1 * mc2) - (mc1 + mc2) + 1 'C 'C Calculate Pearson correlation 'C xmean = 0# ymean = 0# xss = 0# yss = 0# 'C 'C means 'C For i = 1 To mc1 xmean = xmean + i * FMARG(1, i) / total Next i For i = 1 To mc2 ymean = ymean + i * FMARG(2, i) / total Next 'C 'C sums of squared deviations 'C For i = 1 To mc1 xss = xss + FMARG(1, i) * (i - xmean) ^ 2 Next For i = 1 To mc2 yss = yss + FMARG(2, i) * (i - ymean) ^ 2 Next 'C 'C standard deviations 'C xsd = (xss / total) ^ 0.5 ysd = (yss / total) ^ 0.5 'C SUM(x - xbar)(y - ybar) 'C r = ----------------------- 'C ssx ssy scp = 0# For i = 1 To NCELLS scp = scp + FO(i) * ((IVEC(i, 1)) - xmean) * ((IVEC(i, 2)) - ymean) / total scp1 = scp1 + FO(i) * ((IVEC(i, 1)) - xmean) * ((IVEC(i, 1)) - xmean) / total Next corr = scp / (xsd * ysd) prcorr = corr 'C 'C Calculate latent class locations/quadrature points 'C 'C 'C clinc = 10./(LC - 1) 'C b(1) = -5.0 'C do 75 i = 1, lc - 1 'C b(i+1) = -5.0 + (i) * clinc 'C 75 next clinc = 2# * tlimit / (LC - 1) B(1) = -tlimit For i = 1 To LC - 1 B(i + 1) = -tlimit + (i) * clinc Next 'C 'C Calculate normal distribution, if applicable 'C If (IOP(11) = 1) Then Pi = 3.141592653589 const3 = 1# / ((2# * Pi) ^ 0.5) RTOT = 0# For i = 1 To LC r(i) = const3 * Exp(-0.5 * B(i) ^ 2) RTOT = RTOT + r(i) 'C write (2, *) i, b(i), r(i) Next i 'C 'C rescale ordinates to sum to 1.0 'C For i = 1 To LC r(i) = r(i) / RTOT Next i End If 'C 'C Calculate thresholds, if 2-step estimation 'C If (IOP(7) = 1) Then For i = 1 To mc1 - 1 t(1, i + 1) = PPND16(CFMARG(1, i) / total, ier) Next i For i = 1 To mc2 - 1 t(2, i + 1) = PPND16(CFMARG(2, i) / total, ier) Next i End If 'C 'C Calculate default start values 'C If (IOP(1) <> 1) Then 'C 'C Latent class prevalences, if applicable 'C 'C Rho parameter If (corr >= 0#) Then nflag = 0 If (corr > 0.000001) Then x(NEP + 1) = corr Else x(NEP + 1) = 0.002 End If Else nflag = 1 If (corr < -0.000001) Then x(NEP + 1) = -corr Else x(NEP + 1) = 0.002 End If End If 'C 'C Thresholds, if simultaneous estimation 'C INDX = NEP + 1 If (IOP(7) <> 1) Then 'C differing thresholds 'C 'C note: in polycorr, iop(6) is never 1 'C If (IOP(6) <> 1) Then For i = 1 To mc1 - 1 INDX = INDX + 1 x(INDX) = PPND16(CFMARG(1, i) / total, ier) Next i For i = 1 To mc2 - 1 INDX = INDX + 1 x(INDX) = PPND16(CFMARG(2, i) / total, ier) Next i 'C equal thresholds (not applicable to polycorr) Else For i = 1 To MC - 1 INDX = INDX + 1 x(INDX) = PPND16(CFMARG(3, i) / total, ier) Next i End If End If End If 'C 'C Number of estimated parameters 'C 'C NV = NT*(MC - 1) + (2 * LC) + NDIS 'C if (iop(11) <> 0) nv = nv - 2 * lc 'C IF (IOP(6) = 1) NV = NV + (NR - 1) 'C NEST = NV 'C IF (IOP(7) = 1) THEN 'C NV = NV - 2 * (MC - 1) 'C END IF NV = NEP + 1 + (mc1 - 1) + (mc2 - 1) If (IOP(6) = 1) Then NV = NV - (mc2 - 1) NEST = NV 'C 'C adjust if 2-step estimation 'C If (IOP(7) = 1) Then NV = NEP + 1 End If 'C 'C ?? some question about above line 'C 'C 'C User-supplied start values 'C ' 今回は使わない ' if (iop(1) = 1) then ' CALL STARTV (0) 'C 'C if 2-stage estimation, override any user-supplied thresholds 'C ' if (iop(7) = 1) then ' indx = nep + 1 ' for 120 i = 1 to mc1 - 1 ' indx = indx + 1 ' x(indx) = ppnd16(cfmarg(1, i) / total, ier) ' 120 next ' for 121 i = 1 to mc2 - 1 ' indx = indx + 1 ' x(indx) = ppnd16(cfmarg(2, i) / total, ier) ' 121 next ' end if ' end if End Sub Private Sub DATAIN() 'C----------------------------------------------------------------------- 'C SUBROUTINE DATAIN -- INPUT DATA 'C 'C Revised: Change(s) 'C -------- ---------------------------------------------------------- 'C 04/24/92 Made separate subroutine 'C----------------------------------------------------------------------- 'C 'C initialize meta-cell frequencies to 0. 'C ' for i = 1 to nmeta ' go1(i) = 0. ' ge(i) = 0. ' next i total = 0# 'データをここで作成しておく If mc1 = 0 Then mc1 = 2 mc2 = 3 FO(1) = 15 FO(2) = 30 FO(3) = 25 FO(4) = 20 FO(5) = 10 FO(6) = 6 End If NCELLS = mc1 * mc2 'データを直列で「読み込み」 あとで修正 If (IOP(2) = 1) Then GoTo 10 'index frequencies = no =0 固定 ' READ (5, *) (FO(I), I = 1, NCELLS) '//// 'C 'C if there are metacells, read in a pattern array 'C 'これは今回サポートしない ' if (nmeta > 0) then ' read (5, *) (mask2(i), i = 1, ncells) ' end if GoTo 30 10 'READ (1, *) NCELLS ' DO 20 I = 1, NCELLS ' READ (5, *) (IVEC(I, J), J = 1, NR), FO(I) ' 20 CONTINUE 30 total = 0# For i = 1 To NCELLS total = total + FO(i) 'C 'C make sure meta-cell indicator within range 'C If (mask2(i) < 0) Then mask2(i) = 0 If (mask2(i) > nmeta) Then mask2(i) = 0 'C 'C if cell belongs to a meta-cell, update the meta-cell observed 'C frequency 'C If (mask2(i) > 0) Then mm = mask2(i) go1(mm) = go1(mm) + FO(i) End If Next i End Sub Function phi(z) As Double 'C--6/23/00-------------------------------------------------------------- 'C Call appropriate function to calculate normal cdf value 'C 'C If IOP(4) = 0 use alnorm 'C = 1 use normp 'C----------------------------------------------------------------------- If (IOP(4) = 1) Then 'Call normp(z, Prob, q, PDF) 'phi = Prob phi = Application.NormSDist(z) Else If (IOP(4) = 2) Then Call NPROB(z, Prob, q, PDF) phi = Prob Else phi = alnorm(z, False) End If End If End Function Sub NPROB(z, P, q, PDF) ' IMPLICIT DOUBLE PRECISION (A-H,O-Z) 'C 'C P, Q = PROBABILITIES TO THE LEFT AND RIGHT OF Z 'C FOR THE STANDARD NORMAL DISTRIBUTION. 'C PDF = THE PROBABILITY DENSITY FUNCTION 'C 'C REFERENCE: ADAMS,A.G. AREAS UNDER THE NORMAL CURVE, 'C ALGORITHM 39, COMPUTER J., VOL. 12, 197-8, 1969. 'C 'C LATEST REVISION - 23 JANUARY 1981 'C 'C******************************************************************** 'C A0 = 0.5 a1 = 0.398942280444 a2 = 0.399903438504 a3 = 5.75885480458 A4 = 29.8213557808 A5 = 2.62433121679 A6 = 48.6959930692 A7 = 5.92885724438 B0 = 0.398942280385 b1 = 0.000000038052 b2 = 1.00000615302 B3 = 0.000398064794 B4 = 1.98615381364 B5 = 0.151679116635 B6 = 5.29330324926 B7 = 4.8385912808 B8 = 15.1508972451 B9 = 0.742380924027 B10 = 30.789933034 B11 = 3.99019417011 'C ZABS = Abs(z) If (ZABS > 12.7) Then GoTo 20 y = A0 * z * z PDF = Exp(-y) * B0 If (ZABS > 1.28) Then GoTo 10 'C 'C Z BETWEEN -1.28 AND +1.28 'C ' Q = A0 - ZABS * (A1 - A2 * Y / (Y + A3 - A4 / (Y + A5 + A6 / (Y + A7)))) q = A4 / (y + A5 + A6 / (y + A7)) q = A0 - ZABS * (a1 - a2 * y / (y + a3 - q)) If (z < 0#) Then GoTo 30 P = 1# - q Exit Sub 'C 'C ZABS BETWEEN 1.28 AND 12.7 'C 10 q = ZABS - B5 + B6 / (ZABS + B7 - B8 / (ZABS + B9 + B10 / (ZABS + B11))) q = PDF / (ZABS - b1 + b2 / (ZABS + B3 + B4 / q)) If (z < 0#) Then GoTo 30 P = 1# - q Exit Sub 'C 'C Z FAR OUT IN TAIL 'C 20 q = 0# PDF = 0# If (z < 0#) Then GoTo 30 P = 1# Exit Sub 'C 'C NEGATIVE Z, INTERCHANGE P AND Q 'C 30 P = q q = 1# - P Exit Sub End Sub Function PPND16(P, ifault) As Double 'C 'C ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3 'C 'C Produces the normal deviate Z corresponding to a given lower 'C tail area of P; Z is accurate to about 1 part in 10**16. 'C 'C The hash sums below are the sums of the mantissas of the 'C coefficients. They are included for use in checking 'C transcription. 'C Dim r zero = 0# one = 1# half = 0.5 SPLIT1 = 0.425 SPLIT2 = 5# CONST1 = 0.180625 CONST2 = 1.6 'C 'C Coefficients for P close to 0.5 'C A0 = 3.38713287279637 a1 = 133.141667891784 a2 = 1971.59095030655 a3 = 13731.6937655095 A4 = 45921.9539315499 A5 = 67265.7709270087 A6 = 33430.5755835881 A7 = 2509.08092873012 b1 = 42.3133307016009 b2 = 687.187007492058 B3 = 5394.19602142475 B4 = 21213.7943015866 B5 = 39307.8958000927 B6 = 28729.0857357219 B7 = 5226.49527885285 'C HASH SUM AB 55.88319 28806 14901 4439 'C 'C Coefficients for P not close to 0, 0.5 or 1. 'C C0 = 1.42343711074968 c1 = 4.63033784615655 c2 = 5.76949722146069 c3 = 3.6478483247632 c4 = 1.27045825245237 c5 = 0.241780725177451 c6 = 2.27238449892692E-02 C7 = 7.74545014278341E-04 d1 = 2.05319162663776 d2 = 1.6763848301838 d3 = 0.6897673349851 d4 = 0.14810397642748 d5 = 1.51986665636165E-02 D6 = 5.47593808499535E-04 D7 = 1.05075007164442E-09 'C HASH SUM CD 49.33206 50330 16102 89036 'C 'C Coefficients for P near 0 or 1. 'C E0 = 6.6579046435011 E1 = 5.46378491116411 E2 = 1.78482653991729 E3 = 0.296560571828505 E4 = 2.65321895265761E-02 E5 = 1.24266094738808E-03 E6 = 2.71155556874349E-05 E7 = 2.01033439929229E-07 F1 = 0.599832206555888 F2 = 0.136929880922736 F3 = 1.48753612908506E-02 F4 = 7.86869131145613E-04 F5 = 1.84631831751005E-05 F6 = 1.42151175831645E-07 F7 = 2.04426310338994E-15 'C HASH SUM EF 47.52583 31754 92896 71629 'C ifault = 0 q = P - half If (Abs(q) <= SPLIT1) Then r = CONST1 - q * q PPND16 = q * (((((((A7 * r + A6) * r + A5) * r + A4) * r + a3) * r + a2) * r + a1) * r + A0) PPND16 = PPND16 / (((((((B7 * r + B6) * r + B5) * r + B4) * r + B3) * r + b2) * r + b1) * r + one) ' PPND16 = Q * (((((((A7 * R + A6) * R + A5) * R + A4) * R + A3) * R + A2) * R + A1) * R + A0) / (((((((B7 * R + B6) * R + B5) * R + B4) * R + B3) * R + B2) * R + B1) * R + ONE) Exit Function Else If (q < zero) Then r = P Else r = one - P End If If (r <= zero) Then ifault = 1 PPND16 = zero Return End If r = Sqr(-Log(r)) If (r <= SPLIT2) Then r = r - CONST2 PPND16 = (((((((C7 * r + c6) * r + c5) * r + c4) * r + c3) * r + c2) * r + c1) * r + C0) PPND16 = PPND16 / (((((((D7 * r + D6) * r + d5) * r + d4) * r + d3) * r + d2) * r + d1) * r + one) Else r = r - SPLIT2 PPND16 = (((((((E7 * r + E6) * r + E5) * r + E4) * r + E3) * r + E2) * r + E1) * r + E0) PPND16 = PPND16 / (((((((F7 * r + F6) * r + F5) * r + F4) * r + F3) * r + F2) * r + F1) * r + one) End If If (q < zero) Then PPND16 = -PPND16 Exit Function End If End Function Function alnorm(x As Double, upper As Boolean) As Double 'false の場合, =NORMSDIST(x) と同じ 'C 'C Algorithm AS66 Applied Statistics (1973) vol22 no.3 'C 'C Evaluates the tail area of the standardised normal curve 'C from x to infinity if upper is .true. or 'C from minus infinity to x if upper is .false. 'C Dim zero As Double, one As Double, half As Double Dim con As Double, z As Double, y As Double Dim P As Double, q As Double, r As Double Dim a1 As Double, a2 As Double, a3 As Double Dim b1 As Double, b2 As Double Dim c1 As Double, c2 As Double, c3 As Double, c4 As Double, c5 As Double, c6 As Double Dim d1 As Double, d2 As Double, d3 As Double, d4 As Double, d5 As Double Dim up As Boolean 'C*** machine dependent constants Dim ltone As Double, utzero As Double zero = 0# one = 1# half = 0.5 ltone = 7# utzero = 18.66 con = 1.28 P = 0.398942280444 q = 0.39990348504 r = 0.398942280385 a1 = 5.75885480458 a2 = 2.62433121679 a3 = 5.92885724438 b1 = -29.8213557807 b2 = 48.6959930692 c1 = -0.000000038052 c2 = 0.000398064794 c3 = -0.151679116635 c4 = 4.8385912808 c5 = 0.742380924027 c6 = 3.99019417011 d1 = 1.00000615302 d2 = 1.98615381364 d3 = 5.29330324926 d4 = -15.1508972451 d5 = 30.789933034 'C up = upper z = x If (z >= zero) Then GoTo 10 up = Not up z = -z 10 If (z <= ltone Or up And z <= utzero) Then GoTo 20 alnorm = zero GoTo 40 20 y = half * z * z If (z > con) Then GoTo 30 'C alnorm = (P - q * y / (y + a1 + b1 / (y + a2 + b2 / (y + a3)))) alnorm = z * alnorm alnorm = half - alnorm Debug.Print alnorm GoTo 40 30 alnorm = (z + c3 + d3 / (z + c4 + d4 / (z + c5 + d5 / (z + c6)))) alnorm = r * Exp(-y) / (z + c1 + d1 / (z + c2 + d2 / alnorm)) 40 If (Not up) Then alnorm = one - alnorm End Function