program MAIN
USE NUMERICAL_LIBRARIES
implicit none
!==================================== 声明变量 ====================================!
real*8 dClkStart, dClkFinish
integer i, j, k, p, q, r, iRank
integer nFirst, nSecond, nSimu, nMaxRL
parameter (nFirst=1000000, nSecond=100, nSimu=10000, nMaxRL=3000)
real*8 Temp1(1,1)
real*8 Temp2(1,1)
real*8 Temp3
! Model
integer d, nCell, d2, d3
parameter (d=3, nCell=36, d2=3, d3=0) ! No. of variables, each variable 3 3 4
integer nCatLevel(d,1)
! Phase I
integer c1, c2, c3
real*8 Thrld(5,d), Temp4(nCell,1)
real*8 P0(nCell,1), Table(nCell,1), P1(nCell,1)
real*8 X(nCell,d+d2+d3), Y(4,d), W(4,d), M(4,d), Y1(4,d)
real*8 Sigma0(nCell,nCell)
real*8 Rdm(d,1), One(nCell,1), beta0(d+d2,1), lbeta1(d+d2,1), lbeta2(d+d2,d+d2), KP(nCell,1)
real*8 First(d, nFirst), Second(d, nSecond)
real*8 Cov0(d,d), Cov1(d,d), RSIG0(d,d), RSIG1(d,d)
! Phase II
real*8 lambda, Size
real*8 n(nCell,1), z(nCell,1)
integer nIndex
real*8 Limit, MaxA, LimitL, LimitR
real*8 dRL(nSimu,1)
real*8 dARL
real*8 dStdRL
30 FORMAT(8F10.3)
31 FORMAT(8F10.0)
open(UNIT=10, file='CAL-020508-02.txt')
call cpu_time(dClkStart)
!==================================== 准备阶段 ====================================!
data nCatLevel /3, 3, 4/
data P0 /nCell*1.0/
data Thrld /-1000, -0.7, 1, 1000, 0, -1000, -1, 0.5, 1000, 0, -1000, -1.2, 0.3, 2, 1000/ !阈值
data Cov0 /1.0, 0.2, 0.5, 0.2, 1.0, 0.8, 0.5, 0.8, 1.0/
call RNSET(700)
call DCHFAC(d, Cov0, d, 100*DMACH(4), iRank, RSIG0, d)
do i = 1,nFirst
call DRNMVN(1, d, RSIG0, d, Rdm, 1)
First(:,i) = Rdm(:,1)
end do
! 利用阈值，将多元正态分布数据转化为多元有序分类数据（列联表数据）
k = 0
do c1 = 1,nCatLevel(1,1)
do c2 = 1,nCatLevel(2,1)
do c3 = 1,nCatLevel(3,1)
k = k+1
Table(k,1) = dble(count(First(1,:)>Thrld(c1,1) .and. First(1,:)<=Thrld(c1+1,1) .and. &
First(2,:)>Thrld(c2,2) .and. First(2,:)<=Thrld(c2+1,2) .and. &
First(3,:)>Thrld(c3,3) .and. First(3,:)<=Thrld(c3+1,3)))
enddo
enddo
enddo
!write(10,*) Table
!write(10,*)
P0 = Table/dble(nFirst) ! 列联表概率估计
write(10,*) "original cell"
write(10,*) P0
!write(10,*) sum(P0)
write(10,*)
! 计算每个分类变量的边际概率
do i = 1,d
do j = 1,nCatLevel(i,1)
Y(j,i) = dble(count(First(i,:)>Thrld(j,i) .and. First(i,:)<=Thrld(j+1,i)))
end do
end do
Y = Y/dble(nFirst)
write(10,*) "original marginal"
write(10,*) Y
write(10,*)
! 计算每个变量的累积概率
do i = 1,d
do j = 2,nCatLevel(i,1)
Y(j,i) = Y(j-1,i) + Y(j,i)
end do
end do
!write(10,*) Y
!write(10,*)
! score
W(1,:) = Y(1,:)-1
do i = 2,4
W(i,:) = Y(i-1,:)+Y(i,:)-1
end do
write(10,*) "original score"
write(10,*) W
write(10,*)
! 设计主效应相关的矩阵
do i = 1,d
M = W
do j = 1,d
if (j /= i) then
M(:,j) = 1
end if
end do
call KRONECKER(M, nCatLevel, X(:,i)) !Kronecker计算
end do
! 设计交互效应相关的矩阵
do i = 1,d2
M = W
M(:,d+1-i) = 1
call KRONECKER(M, nCatLevel, X(:,i+d)) !Kronecker计算
end do
Sigma0 = 0
call DIAGMAT(P0, Sigma0)
Sigma0 = Sigma0 - (P0 .x. .t.P0)
! paramters
lambda = 0.2
Size = dble(nSecond)
LimitL = 2
LimitR = 2.6
!==================================== 计算统计量的值（nSimu=10000个） ====================================!
call RNSET(12545)
dRL = 0
nIndex = 0
MaxA = 0
z = P0*Size
do i = 1,nSimu
Second = 0
do k =1,nSecond
call DRNMVN(1, d, RSIG0, d, Rdm, 1) !仿真出样本
Second(:,k) = Rdm(:,1)
enddo
r = 0 !分割成分类数据
do c1 = 1,nCatLevel(1,1)
do c2 = 1,nCatLevel(2,1)
do c3 = 1,nCatLevel(3,1)
r = r+1
n(r,1) = dble(count(Second(1,:)>Thrld(c1,1) .and. Second(1,:)<=Thrld(c1+1,1) .and. &
Second(2,:)>Thrld(c2,2) .and. Second(2,:)<=Thrld(c2+1,2) .and. &
Second(3,:)>Thrld(c3,3) .and. Second(3,:)<=Thrld(c3+1,3)))
enddo
enddo
enddo
Temp1 = .t.(n-P0*Size) .x. X .x. .i.(.t.X .x. Sigma0 .x. X) .x. .t.X .x. (n-P0*Size) !计算统计量的值
MaxA = Temp1(1,1) / Size
dRL(nIndex,1) = MaxA !储存统计量的值
write(*,*) i,j,MaxA
end do
call cpu_time(dClkFinish)
write(*,*) "Time is", dClkFinish-dClkStart, "seconds"
stop
end
subroutine DIAGMAT(Elm, Mat)
implicit none
integer nCell
parameter(nCell = 36)
real*8 Elm(nCell,1)
real*8 Mat(nCell,nCell)
integer i
do i = 1,nCell
Mat(i,i) = Elm(i,1)
end do
return
end
subroutine KRONECKER(M, nDim, MK) !!Kronecker计算
implicit none
integer i, j, p, q
real*8 M(4,3)
integer nDim(3,1)
real*8 MK(36,1)
real*8,allocatable :: MK1(:,:)
real*8,allocatable :: M1(:,:)
real*8,allocatable :: M2(:,:)
real*8,allocatable :: M3(:,:)
30 FORMAT(6F10.0)
if (.not. allocated(M1)) then
allocate(M1(nDim(1,1),1))
end if
if (.not. allocated(M2)) then
allocate(M2(nDim(2,1),1))
end if
if (.not. allocated(M3)) then
allocate(M3(nDim(3,1),1))
end if
if (.not. allocated(MK1)) then
allocate(MK1(nDim(1,1)*nDim(2,1),1))
end if
M1(:,1) = M(1:nDim(1,1),1)
M2(:,1) = M(1:nDim(2,1),2)
M3(:,1) = M(1:nDim(3,1),3)
do i = 1,nDim(1,1)
do j = 1,nDim(2,1)
MK1((i-1)*nDim(2,1)+j,1) = M1(i,1)*M2(j,1)
end do
end do
do i = 1,nDim(1,1)*nDim(2,1)
do j = 1,nDim(3,1)
MK((i-1)*nDim(3,1)+j,1) = MK1(i,1)*M3(j,1)
end do
end do
return
end