(*  Mathematica Routines for Statistical Models for Behavioral Observations
            David Rogosa and Ghassan Ghandour
                  November 1991
            Copyright (c)  1991

This file contains Mathematica Expressions Corresponding to most
of the equations in Rogosa and Ghandour 1991, Statistical Models for
Behavioral Observations *)

(* Some Notational Compromises
   Equation numbering:  Equation 1.2' p.174 is indicated as
   eq1p2prime , and so forth
   Conditioning: $ stands in for the | character in conditional 
   expections and variances (e.g. eq1p3 below)
   A rendition of the Glossary pp. 158-9, adapted for the notation
   of these routines, is appended to the end of this file in comment
   form.
  *)

(* Useful Constants Vectors
To keep in touch with your settings for parameter values,
we provide 4 vectors of commonly used values that can be called
up during calculations-- indConst, popConst, recConst, timeConst
*)
indConst := {"alpha = ",alpha, "inc = ",inc,  "phi = ",phi,
             "Iprocess = ",Iprocess, "mu1 = ",mu1,
             "tau2 = ",tau2, "T = ",T}

popConst := {"mualpha = ",mualpha, "sig2alpha = ",sig2alpha,
             "muinc = ",muinc, "sig2inc = ",sig2inc,
             "corrincalpha = ",corrincalpha, "Iprocess = ",Iprocess,
             "muphi = ",muphi,"sig2phi = ",sig2phi,
             "corrincphi = ",corrincphi, "mumu1 = ",mumu1,
             "sig2mu1 = ",sig2mu1}

recConst := {"R = ",R, "kap = ",kap, "gam = ",gam, "om = ",om,
             "psi = ", psi,
             "minkap = ",minkap, "rangekap = ",rangekap,
             "mingam = ",mingam, "rangegam = ",rangegam,
             "corrkapgam = ",corrkapgam, "muom = ",muom,
             "sig2om = ",sig2om, "mupsi = ",mupsi, "mupsif = ",mupsif}

timeConst := {"occ = ",occ, "T = ",T, "s = ",s, "s2 = ",s2, "al = ",al,
              "tau2 = ",tau2, "deltaInMinutes = ",delta*60, 
              "segs = ",T/delta}


eq1p2prime ::= vpT 

vpT := (inc^3 sig2w)/T + (1/6 + (inc^4 sig2w^2)/2 - 
                                (inc^3 mu3w)/3)/(T^2)

(* For some calculations use a function form for vpT *)

funcvpT[T_, inc_, sig2w_, mu3w_] := (inc^3 sig2w)/T + 
                                    (1/6 + (inc^4 sig2w^2)/2 - 
                                    (inc^3 mu3w)/3)/(T^2)

(* For gamma point process sig2w is  1/(alpha inc^2)
   mu3w  is  2/(alpha^2 inc^3)   *)
(* For alternating Poisson sig2w is (1 - 2*phi + 2*phi^2)/(inc^2)
   mu3w is (2 - 6*phi + 6*phi^2)/(inc^3)  *)


(* Recorder parameters and derived quantities, single recorder *)
om := 1 - kap + gam
psi := kap*(1 - kap) + gam

(* Moments for a fixed recorder *)

eq1p3a ::= EMprDivT$pr 

EMprDivT$pr := inc om


eq1p3b ::= VMprDivT$pr

VMprDivT$pr := (inc psi)/T + om^2 vpT

(* Recorder parameters and derived quantities for a population of
recorders, assuming uniform distributions for kap and gam *)
(* Input Recorder population parameters *)
minkap = .
rangekap = .
mingam = .
rangegam = .
corrkapgam = .

(* Standard Omega/Psi Moments from kappa gamma Uniform *)

Clear[muunf]
Clear[varunf]
Clear[meanpsi]
Clear[meanom]
Clear[varomg]
      
(* Uniform Moments *)
muunf[min_, range_] := (range + 2 min) / 2
varunf[range_] := range^2 / 12

(* Theoretical Moments  Page 176  Applied to Uniform Distributions *)
meanpsi[minkap_, rangekap_, mingam_, rangegam_] := 
                 muunf[minkap, rangekap] - varunf[rangekap] - 
                 muunf[minkap, rangekap]^2 + muunf[mingam, rangegam]

meanom[minkap_, rangekap_, mingam_, rangegam_] := 
                 1 -  muunf[minkap, rangekap] + muunf[mingam, rangegam]

varom[rangekap_, rangegam_, corrkapgam_] := 
                 varunf[rangekap] + varunf[rangegam] -
                 2 corrkapgam Sqrt[varunf[rangekap] varunf[rangegam]]


muom := meanom[minkap, rangekap, mingam, rangegam]
sig2om := varom[rangekap, rangegam,corrkapgam]
mupsi := meanpsi[minkap, rangekap, mingam, rangegam]

(* For fixed recorder we use muspif to represent degenerate 
   kappa gamma distribution *)
mupsif := meanpsi[minkap, rangekap, mingam, rangegam] +  varunf[rangekap] 

(* Equation 1.4 *)

eq1p4a ::= EMprDivT$p

EMprDivT$p := inc muom

eq1p4b ::= VMprDivT$p

VMprDivT$p := (inc mupsi/T) + (muom^2 + sig2om)*vpT + 
              inc^2 sig2om

(* Equation 1.6 *)
eq1p6a ::= EMbarpDivT$p

EMbarpDivT$p := inc muom

eq1p6b ::= VMbarpDivT$p

VMbarpDivT$p := (muom^2 vpT) + (1/R)((inc mupsi/T) + 
                sig2om*(inc^2 + vpT))


eq1p7b ::= VNpoDivT$p

VNpoDivT$p := tau2 + vpT

(* Defining Observation Schedule S and Allocation Matrix A *)
Clear[kp]
Clear[alKeep]
Clear[alDiag]
Clear[s]
Clear[one]
Clear[vpos]
Clear[saatst]
Clear[sst]
Clear[svt]
Clear[svt1]
Clear[svt1st]
Clear[Tplus]

(* Note--- dimension of s, i.e. number of occasions O = occ needs to be
specified before computation *)
occ = 1

(* Allocation matrix for random recorder rosters resampled each Occ *)
ADiag := IdentityMatrix[Round[occ]]

(* AKeep is allocation matrix for keeping same recorder roster 
see pages 178, 180 *)
kp[ik_,jk_] := If[jk == 1, 1, 0, 0]
AKeep := Table[kp[iak, jak], {iak, Round[occ]}, {jak, Round[occ]}]

(* Initial Setting--Allocation Matrix resampled recorder  *)
al := ADiag

(* Observation schedule vectors *)
(* for equal observation intervals set s := Table[T, {is, Round[occ]}] *)

s = {T}
s2 := Table[s[[iss]]^2, {iss, Round[occ]}]


(* Vpo(s) vector *)
vpos := Table[funcvpT[s[[iv]], inc, sig2w, mu3w],{iv, Round[occ]}]

(* Matrix products *)
saatst := s.al.Transpose[al].s
sst := s.s

(* Tp+ *)
one := Table[i1/i1, {i1, Round[occ]}]
Tplus := s.one

(* eq1p8b = VNpplusDivTplus$p is obtained from eq1p11b with single
perfect recorder-- minkap = rangekap = mingam = rangegam = 0 --
and R = 1 *)
(* eq1p9b =  VMpplusrDivTplus$p is obtained from eq1p11b with single
recorder, R = 1 *)

(* Equation 1.11: Multiple Recorders, Multiple Occasions *)

eq1p11a ::= EMbarpplusDivTplus$p

EMbarpplusDivTplus$p := inc muom

eq1p11b ::= VMbarpplusDivTplus$p 

VMbarpplusDivTplus$p := (inc^2*sig2om*saatst)/(R*Tplus^2) +
                        (inc*mupsi)/(R*Tplus) +
                        (R muom^2 + sig2om) * 
                        (sst*tau2 + s2.vpos)/(R*Tplus^2)

(* Section 1.3-- Examples of Renewal Processes *)

eq1p12 ::= vpTpoisson

vpTpoisson := funcvpT[T, inc, 1/inc^2, 2/inc^3] 

eq1p13 ::= vpTgamma

vpTgamma := funcvpT[T, inc, 1/(alpha inc^2), 2/(alpha^2 inc^3)] 

eq1p14 ::= vpTaltP

vpTaltP := funcvpT[T, inc, (1 - 2*phi + 2*phi^2)/(inc^2) , 
           (2 - 6*phi + 6*phi^2)/(inc^3)]

(* SECTION 2 *)
(* input for parameter distributions over population of individuals *)
muinc = .
sig2inc = .
muphi = .
sig2phi = .
corrincphi = .
covincphi := corrincphi Sqrt[sig2phi sig2inc]
mualpha = .
sig2alpha = .
corrincalpha = .
covincalpha := corrincalpha Sqrt[sig2alpha sig2inc]


eq2p4 ::= relNDivT

relNDivT := sig2inc/(sig2inc + EpvpT)

eq2p6 ::= EpvpTgamma

EpvpTgamma := (-mualpha^2 + mualpha^4 - 6*covincalpha*mualpha^2*T + 
              6*mualpha^3*muinc*T - 3*sig2alpha + 
              6*mualpha*muinc*T*sig2alpha)/(6*mualpha^4*T^2)

eq2p7 ::= EpvpTaltP

EpvpTaltP := (2*muphi^2 - 4*muphi^3 + 2*muphi^4 - 2*covincphi*T + 
             muinc*T +  4*covincphi*muphi*T - 2*muinc*muphi*T + 
             2*muinc*muphi^2*T + 2*sig2phi - 12*muphi*sig2phi + 
             12*muphi^2*sig2phi + 2*muinc*T*sig2phi)/T^2

(* In funcEpvpT declare Iprocess = 1 for gamma; 
                                   0 for Alternating Poisson *)

funcEpvpT[T_, muinc_, sig2inc_, muphi_, sig2phi_, covincphi_, mualpha_, 
          sig2alpha_, covincalpha_, Iprocess_] :=
          Iprocess*(-mualpha^2 + mualpha^4 - 6*covincalpha*mualpha^2*T + 
                   6*mualpha^3*muinc*T - 3*sig2alpha + 
                   6*mualpha*muinc*T*sig2alpha)/(6*mualpha^4*T^2) +
          (1 - Iprocess)*(2*muphi^2 - 4*muphi^3 + 2*muphi^4 - 
                   2*covincphi*T + muinc*T +  4*covincphi*muphi*T - 
                   2*muinc*muphi*T + 2*muinc*muphi^2*T + 2*sig2phi - 
                   12*muphi*sig2phi + 12*muphi^2*sig2phi + 
                   2*muinc*T*sig2phi)/T^2

eq2p8 := sig2inc/((muinc/muT)*(1 + (sig2T/muT^2) - covincT/(muinc muT)) +
                                                              sig2inc) 

eq2p9 ::= relMprDivT

relMprDivT := sig2inc/(sig2inc + (muinc mupsi)/(T muom^2) +
              EpvpT + (sig2om/muom^2)*(EpvpT +
              sig2inc + muinc^2))

(* Eq2.10b is 2.12b for occ = 1, R = 1  *)

(* Moments over the population for Multiple Occasions *)

Epvpos := Table[funcEpvpT[s[[iv]], muinc, sig2inc, muphi, sig2phi,
                          covincphi, mualpha, sig2alpha, covincalpha,
                          Iprocess],{iv, Round[occ]}]


eq2p12 ::= EpVMbarpplusDivTplus

EpVMbarpplusDivTplus :=  (muinc*mupsi)/(R*Tplus) +
                         (muinc^2 + sig2inc)*sig2om*saatst/
                         (R*Tplus^2) + (R muom^2 + sig2om) * 
                         (sst*tau2 + s2.Epvpos)/(R*Tplus^2)

(* Tables 3 through 5-- Spearman-Brown Illustrations  *)

actualReliability := sig2inc/(sig2inc + eq2p12/muom^2)
 
(* Compute seedReliability for k = 1 for use in Spearman-Brown *)

SBReliability := k*seedReliability/((k - 1)*seedReliability + 1)

(* Note in Table 5, bottom leftmost entry should be .312, not .296  *)

eq3p3 ::= rhoMM'

rhoMM' := (sig2inc + EpvpT)/(sig2inc + (muinc mupsi)/(T muom^2) +
          EpvpT + (sig2om/muom^2)*(EpvpT + sig2inc + muinc^2))

eq3p8 ::= prVC

prVC := {prVCp, prVCr, prVCpr, prVCeps} 

prVCp     := muom^2 sig2inc + muom^2 EpvpT
prVCr     := sig2om muinc^2
prVCpr    := sig2inc sig2om + sig2om EpvpT
prVCeps   := muinc mupsi/T


(* pxoxr G theory quantities *)
eq3p11 ::= porVC

porVC := {porVCp, porVCo, porVCr, porVCpo, porVCpr, porVCro,
                  porVCpor, porVCeps}

porVCp := muom^2 sig2inc
porVCo = 0
porVCr := sig2om muinc^2
porVCpo := muom^2*(EpvpT + tau2)
porVCpr := sig2inc sig2om
porVCro = 0
porVCpor := sig2om*(EpvpT + tau2)
porVCeps := muinc mupsi/T


eq3p12 ::= porVCabs

porVCabs := porVCo + porVCr + porVCpo + porVCpr + porVCro +
            (porVCpor + porVCeps)

eq3p13 ::= porVCrel

porVCrel := porVCpo + porVCpr + (porVCpor + porVCeps)

eq3p15 ::= Varpdel$r

Varpdel$r := sig2inc*(1 - muom)^2 + (muinc mupsif)/T + muom^2 *
             (EpvpT + tau2) 

(* D study projections-- Multiple Occasions and Recorders *)

eq3p16 ::= porVCabsOR

porVCabsOR :=  porVCo/occ + porVCr/R + porVCpo/occ + porVCpr/R + 
               (porVCro + (porVCpor + porVCeps))/(occ R)

eq3p17 ::= porVCrelOR

porVCrelOR :=  porVCo/occ + porVCpo/occ + porVCpr/R + 
               (porVCro + (porVCpor + porVCeps))/(occ R)


(* Eq3.18a is 2.12 with al = ADiag and Eq3.18b is 2.12 with al = AKeep *)

eq3p19 ::= VarpdelOR$scriptR

VarpdelOR$scriptR := sig2inc*(1 - muom)^2 + (muinc mupsif)/
                     (R T occ) + muom^2*(EpvpT + tau2)/occ

(* Section 4  *)

(* Equation 4.1 *)
eq4p1 ::= EC

EC := T muc inc

(* Equation 4.2 *)
eq4p2 ::= VC

VC := T inc sig2c*(1 - rhocw^2) + T^2*funcvpT[T, inc, sig2w, mu3w]*
      (muc - betacw/inc)^2

(* Equation 4.3 *)

eq4p3 ::= CovCN

CovCN := T^2*(muc - betacw/inc)*funcvpT[T, inc, sig2w, mu3w]

(* Equation 4.4a *)

eq4p4a ::= ECDivN

ECDivN := muc + (betacw*funcvpT[T, inc, sig2w, mu3w])/inc^3

(* Equation 4.4b *)
eq4p4b ::= VCDivN

VCDivN := (sig2c/(T inc))*((1 - rhocw^2) + T*
          funcvpT[T, inc, sig2w, mu3w]*(rhocw^2/(sig2w inc^3)))

(*                     For proportions 
muc = pi
sig2c = pi*(1-pi)
rhocw = rho
betacw = rhocw*Sqrt[sig2c/sig2w]
                                  *)

(* Equation 4.5 is obtained from 4.4b using the proportion settings
   for moments of c,w *)
eq4p5 ::= VZDivN$p

VZDivN$p := -((-1 + pi)*pi*(rho^2 - 2*inc^3*mu3w*rho^2 + 
            6*T*inc^3*sig2w + 3*inc^4*rho^2*sig2w^2))/
            (6*T^2*inc^4*sig2w)

(* Equation 4.6 can be expressed more generally for either an
   Alternating Poisson Process or for a Gamma Process
   Here we give 4.6 for the Gamma Process with alphap = alpha for all p
   Either one can be reduced to Poisson case given in p.228 *)
eq4p6Gam ::= EpvpZDivNG

EpvpZDivNG :=
   (2*covincpi*muinc*rho^2 - 2*alpha^2*covincpi*muinc*rho^2 - 
    4*covincpi*muinc*mupi*rho^2 + 4*alpha^2*covincpi*muinc*mupi*rho^2 - 
    muinc^2*mupi*rho^2 + alpha^2*muinc^2*mupi*rho^2 + muinc^2*mupi^2*rho^2 - 
    alpha^2*muinc^2*mupi^2*rho^2 - 3*mupi*rho^2*sig2inc + 
    3*alpha^2*mupi*rho^2*sig2inc + 3*mupi^2*rho^2*sig2inc - 
    3*alpha^2*mupi^2*rho^2*sig2inc + muinc^2*rho^2*sig2pi - 
    alpha^2*muinc^2*rho^2*sig2pi - 6*alpha*covincpi*muinc^2*T + 
    12*alpha*covincpi*muinc^2*mupi*T + 6*alpha*muinc^3*mupi*T - 
    6*alpha*muinc^3*mupi^2*T + 6*alpha*muinc*mupi*sig2inc*T - 
    6*alpha*muinc*mupi^2*sig2inc*T - 6*alpha*muinc^3*sig2pi*T)/
    (6*alpha*muinc^4*T^2)

(* For Prevalence and Event Duration From Alternating Poisson
muc := phi/inc
sig2c := muc^2
betacw := phi^2 / ((1-phi)^2 + phi^2)
rhocw := Sqrt[betacw]
                                     *)

eq4p7 ::= VLDivT$p

VLDivT$p := ((phi^2*(1 - phi^2/((1 - phi)^2 + phi^2))*T)/inc + 
            (phi/inc - phi^2/(inc*((1 - phi)^2 + phi^2)))^2*
            ((2*(1 - phi)^2*phi^2)/T^2 + (inc*((1 - phi)^2 + 
            phi^2))/T)*T^2)/T^2

(* Section 4.4.3, Psychometric Properties of Empirical Prevalence *)

page231Epvpprev := 
(2*T*muinc^3*muphi^2 - 20*T*muinc^3*muphi^3 + 2*muinc^2*muphi^4 + 
    98*T*muinc^3*muphi^4 - 24*muinc^2*muphi^5 - 304*T*muinc^3*muphi^5 + 
    132*muinc^2*muphi^6 + 656*T*muinc^3*muphi^6 - 440*muinc^2*muphi^7 - 
    1024*T*muinc^3*muphi^7 + 986*muinc^2*muphi^8 + 1168*T*muinc^3*muphi^8 - 
    1552*muinc^2*muphi^9 - 960*T*muinc^3*muphi^9 + 1736*muinc^2*muphi^10 + 
    544*T*muinc^3*muphi^10 - 1360*muinc^2*muphi^11 - 192*T*muinc^3*muphi^11 + 
    712*muinc^2*muphi^12 + 32*T*muinc^3*muphi^12 - 224*muinc^2*muphi^13 + 
    32*muinc^2*muphi^14 + 2*T*muinc*muphi^2*sig2inc - 
    20*T*muinc*muphi^3*sig2inc + 6*muphi^4*sig2inc + 
    98*T*muinc*muphi^4*sig2inc - 72*muphi^5*sig2inc - 
    304*T*muinc*muphi^5*sig2inc + 396*muphi^6*sig2inc + 
    656*T*muinc*muphi^6*sig2inc - 1320*muphi^7*sig2inc - 
    1024*T*muinc*muphi^7*sig2inc + 2958*muphi^8*sig2inc + 
    1168*T*muinc*muphi^8*sig2inc - 4656*muphi^9*sig2inc - 
    960*T*muinc*muphi^9*sig2inc + 5208*muphi^10*sig2inc + 
    544*T*muinc*muphi^10*sig2inc - 4080*muphi^11*sig2inc - 
    192*T*muinc*muphi^11*sig2inc + 2136*muphi^12*sig2inc + 
    32*T*muinc*muphi^12*sig2inc - 672*muphi^13*sig2inc + 
    96*muphi^14*sig2inc - 4*T*corrincphi*muinc^2*muphi*sig2inc^(1/2)*
     sig2phi^(1/2) + 44*T*corrincphi*muinc^2*muphi^2*sig2inc^(1/2)*
     sig2phi^(1/2) - 16*corrincphi*muinc*muphi^3*sig2inc^(1/2)*
     sig2phi^(1/2) - 232*T*corrincphi*muinc^2*muphi^3*sig2inc^(1/2)*
     sig2phi^(1/2) + 208*corrincphi*muinc*muphi^4*sig2inc^(1/2)*
     sig2phi^(1/2) + 768*T*corrincphi*muinc^2*muphi^4*sig2inc^(1/2)*
     sig2phi^(1/2) - 1200*corrincphi*muinc*muphi^5*sig2inc^(1/2)*
     sig2phi^(1/2) - 1760*T*corrincphi*muinc^2*muphi^5*sig2inc^(1/2)*
     sig2phi^(1/2) + 4112*corrincphi*muinc*muphi^6*sig2inc^(1/2)*
     sig2phi^(1/2) + 2912*T*corrincphi*muinc^2*muphi^6*sig2inc^(1/2)*
     sig2phi^(1/2) - 9376*corrincphi*muinc*muphi^7*sig2inc^(1/2)*
     sig2phi^(1/2) - 3520*T*corrincphi*muinc^2*muphi^7*sig2inc^(1/2)*
     sig2phi^(1/2) + 14976*corrincphi*muinc*muphi^8*sig2inc^(1/2)*
     sig2phi^(1/2) + 3072*T*corrincphi*muinc^2*muphi^8*sig2inc^(1/2)*
     sig2phi^(1/2) - 17056*corrincphi*muinc*muphi^9*sig2inc^(1/2)*
     sig2phi^(1/2) - 1856*T*corrincphi*muinc^2*muphi^9*sig2inc^(1/2)*
     sig2phi^(1/2) + 13728*corrincphi*muinc*muphi^10*sig2inc^(1/2)*
     sig2phi^(1/2) + 704*T*corrincphi*muinc^2*muphi^10*sig2inc^(1/2)*
     sig2phi^(1/2) - 7488*corrincphi*muinc*muphi^11*sig2inc^(1/2)*
     sig2phi^(1/2) - 128*T*corrincphi*muinc^2*muphi^11*sig2inc^(1/2)*
     sig2phi^(1/2) + 2496*corrincphi*muinc*muphi^12*sig2inc^(1/2)*
     sig2phi^(1/2) - 384*corrincphi*muinc*muphi^13*sig2inc^(1/2)*
     sig2phi^(1/2) + 2*T*muinc^3*sig2phi - 28*T*muinc^3*muphi*sig2phi + 
    12*muinc^2*muphi^2*sig2phi + 172*T*muinc^3*muphi^2*sig2phi - 
    176*muinc^2*muphi^3*sig2phi - 640*T*muinc^3*muphi^3*sig2phi + 
    1084*muinc^2*muphi^4*sig2phi + 1616*T*muinc^3*muphi^4*sig2phi - 
    3832*muinc^2*muphi^5*sig2phi - 2912*T*muinc^3*muphi^5*sig2phi + 
    8888*muinc^2*muphi^6*sig2phi + 3808*T*muinc^3*muphi^6*sig2phi - 
    14464*muinc^2*muphi^7*sig2phi - 3584*T*muinc^3*muphi^7*sig2phi + 
    16936*muinc^2*muphi^8*sig2phi + 2336*T*muinc^3*muphi^8*sig2phi - 
    14160*muinc^2*muphi^9*sig2phi - 960*T*muinc^3*muphi^9*sig2phi + 
    8112*muinc^2*muphi^10*sig2phi + 192*T*muinc^3*muphi^10*sig2phi - 
    2880*muinc^2*muphi^11*sig2phi + 480*muinc^2*muphi^12*sig2phi)/
  (T^2*muinc^4*(1 - 2*muphi + 2*muphi^2)^4)

relprev := sig2phi/(page231Epvpprev + sig2phi)

(* Section 4.5  Event Duration  *)
(* parameter value *)
mu1 := phi/inc

page232biasDur ::= biasLDivN$p

biasLDivN$p :=
(phi^2*(2*phi^2 - 4*phi^3 + 2*phi^4 + inc*T - 2*inc*phi*T + 
2*inc*phi^2*T))/(inc^3*(1 - 2*phi + 2*phi^2)*T^2)

eq4p8 ::= VLDivN$p

VLDivN$p :=
(phi^2*(2*phi^4 - 4*phi^5 + 2*phi^6 + inc*T - 4*inc*phi*T + 
8*inc*phi^2*T - 8*inc*phi^3*T + 4*inc*phi^4*T))/(inc^4*(1 - 2*phi + 
2*phi^2)^2*T^2)  

page232Epvpdur :=
(T*muinc^3*muphi^2 - 8*T*muinc^3*muphi^3 + 32*T*muinc^3*muphi^4 - 
    80*T*muinc^3*muphi^5 + 2*muinc^2*muphi^6 + 136*T*muinc^3*muphi^6 - 
    12*muinc^2*muphi^7 - 160*T*muinc^3*muphi^7 + 34*muinc^2*muphi^8 + 
    128*T*muinc^3*muphi^8 - 56*muinc^2*muphi^9 - 64*T*muinc^3*muphi^9 + 
    56*muinc^2*muphi^10 + 16*T*muinc^3*muphi^10 - 32*muinc^2*muphi^11 + 
    8*muinc^2*muphi^12 + 6*T*muinc*muphi^2*sig2inc - 
    48*T*muinc*muphi^3*sig2inc + 192*T*muinc*muphi^4*sig2inc - 
    480*T*muinc*muphi^5*sig2inc + 20*muphi^6*sig2inc + 
    816*T*muinc*muphi^6*sig2inc - 120*muphi^7*sig2inc - 
    960*T*muinc*muphi^7*sig2inc + 340*muphi^8*sig2inc + 
    768*T*muinc*muphi^8*sig2inc - 560*muphi^9*sig2inc - 
    384*T*muinc*muphi^9*sig2inc + 560*muphi^10*sig2inc + 
    96*T*muinc*muphi^10*sig2inc - 320*muphi^11*sig2inc + 
    80*muphi^12*sig2inc - 6*T*corrincphi*muinc^2*muphi*sig2inc^(1/2)*
     sig2phi^(1/2) + 48*T*corrincphi*muinc^2*muphi^2*sig2inc^(1/2)*
     sig2phi^(1/2) - 192*T*corrincphi*muinc^2*muphi^3*sig2inc^(1/2)*
     sig2phi^(1/2) + 480*T*corrincphi*muinc^2*muphi^4*sig2inc^(1/2)*
     sig2phi^(1/2) - 48*corrincphi*muinc*muphi^5*sig2inc^(1/2)*
     sig2phi^(1/2) - 816*T*corrincphi*muinc^2*muphi^5*sig2inc^(1/2)*
     sig2phi^(1/2) + 272*corrincphi*muinc*muphi^6*sig2inc^(1/2)*
     sig2phi^(1/2) + 960*T*corrincphi*muinc^2*muphi^6*sig2inc^(1/2)*
     sig2phi^(1/2) - 704*corrincphi*muinc*muphi^7*sig2inc^(1/2)*
     sig2phi^(1/2) - 768*T*corrincphi*muinc^2*muphi^7*sig2inc^(1/2)*
     sig2phi^(1/2) + 1056*corrincphi*muinc*muphi^8*sig2inc^(1/2)*
     sig2phi^(1/2) + 384*T*corrincphi*muinc^2*muphi^8*sig2inc^(1/2)*
     sig2phi^(1/2) - 960*corrincphi*muinc*muphi^9*sig2inc^(1/2)*
     sig2phi^(1/2) - 96*T*corrincphi*muinc^2*muphi^9*sig2inc^(1/2)*
     sig2phi^(1/2) + 512*corrincphi*muinc*muphi^10*sig2inc^(1/2)*
     sig2phi^(1/2) - 128*corrincphi*muinc*muphi^11*sig2inc^(1/2)*
     sig2phi^(1/2) + T*muinc^3*sig2phi - 8*T*muinc^3*muphi*sig2phi + 
    32*T*muinc^3*muphi^2*sig2phi - 80*T*muinc^3*muphi^3*sig2phi + 
    30*muinc^2*muphi^4*sig2phi + 136*T*muinc^3*muphi^4*sig2phi - 
    156*muinc^2*muphi^5*sig2phi - 160*T*muinc^3*muphi^5*sig2phi + 
    344*muinc^2*muphi^6*sig2phi + 128*T*muinc^3*muphi^6*sig2phi - 
    448*muinc^2*muphi^7*sig2phi - 64*T*muinc^3*muphi^7*sig2phi + 
    376*muinc^2*muphi^8*sig2phi + 16*T*muinc^3*muphi^8*sig2phi - 
    192*muinc^2*muphi^9*sig2phi + 48*muinc^2*muphi^10*sig2phi)/
  (T^2*muinc^6*(1 - 2*muphi + 2*muphi^2)^4)

(* Moments of Event Duration for the {Inc, Phi} distribution *)
mumu1 :=  muphi/muinc - (1/muinc^2)*covincphi + (muphi/muinc^3)*sig2inc

sig2mu1 := (muphi/muinc)^2*(sig2phi/muphi^2 + sig2inc/muinc^2 -
                          2 covincphi/(muphi muinc))

reldur := sig2mu1/(page232Epvpdur + sig2mu1)

        (* Section 4.6   Relative Precisions  *)
(* CV for individual p , Section 4.6.1  *)

incCV$p :=  Sqrt[vpT]/inc

prevCV$p := Sqrt[eq4p7]/phi

durCV$p  := Sqrt[eq4p8]/mu1

                          (*  Section 5  *)
(*  define segs = T/delta  *)

eq5p1 ::= VelDivsegs$p

VelDivsegs$p := (delta phi*(1 - phi)/T)* 
                (2/(1 - Exp[-inc delta/(phi*(1 - phi))]) - 1)

(* Entries for Table 8 *)

table8entries := ((-1 + 2/(1 - E^((delta*inc)/((-1 + phi)*phi))))*delta*
                inc^2*(1 - phi)*(1 - 2*phi + 2*phi^2)^2*T)/
                (2*(-1 + phi)^2*phi*(phi^2 - 6*phi^3 + 13*phi^4 - 
                12*phi^5 + 4*phi^6 + inc*T - 4*inc*phi*T + 
                8*inc*phi^2*T - 8*inc*phi^3*T + 4*inc*phi^4*T))

(* Ratio of Instantaneous to Continuous reliabilities, page 237 *)

page237RelRatio := (sig2phi + (2*T*muinc^3*muphi^2 - 
       20*T*muinc^3*muphi^3 + 2*muinc^2*muphi^4 + 
       98*T*muinc^3*muphi^4 - 24*muinc^2*muphi^5 - 304*T*muinc^3*muphi^5 + 
       132*muinc^2*muphi^6 + 656*T*muinc^3*muphi^6 - 440*muinc^2*muphi^7 - 
       1024*T*muinc^3*muphi^7 + 986*muinc^2*muphi^8 + 
       1168*T*muinc^3*muphi^8 - 1552*muinc^2*muphi^9 - 
       960*T*muinc^3*muphi^9 + 1736*muinc^2*muphi^10 + 
       544*T*muinc^3*muphi^10 - 1360*muinc^2*muphi^11 - 
       192*T*muinc^3*muphi^11 + 712*muinc^2*muphi^12 + 
       32*T*muinc^3*muphi^12 - 224*muinc^2*muphi^13 + 32*muinc^2*muphi^14 + 
       2*T*muinc*muphi^2*sig2inc - 20*T*muinc*muphi^3*sig2inc + 
       6*muphi^4*sig2inc + 98*T*muinc*muphi^4*sig2inc - 72*muphi^5*sig2inc - 
       304*T*muinc*muphi^5*sig2inc + 396*muphi^6*sig2inc + 
       656*T*muinc*muphi^6*sig2inc - 1320*muphi^7*sig2inc - 
       1024*T*muinc*muphi^7*sig2inc + 2958*muphi^8*sig2inc + 
       1168*T*muinc*muphi^8*sig2inc - 4656*muphi^9*sig2inc - 
       960*T*muinc*muphi^9*sig2inc + 5208*muphi^10*sig2inc + 
       544*T*muinc*muphi^10*sig2inc - 4080*muphi^11*sig2inc - 
       192*T*muinc*muphi^11*sig2inc + 2136*muphi^12*sig2inc + 
       32*T*muinc*muphi^12*sig2inc - 672*muphi^13*sig2inc + 
       96*muphi^14*sig2inc - 4*T*corrincphi*muinc^2*muphi*sig2inc^(1/2)*
        sig2phi^(1/2) + 44*T*corrincphi*muinc^2*muphi^2*sig2inc^(1/2)*
        sig2phi^(1/2) - 16*corrincphi*muinc*muphi^3*sig2inc^(1/2)*
        sig2phi^(1/2) - 232*T*corrincphi*muinc^2*muphi^3*sig2inc^(1/2)*
        sig2phi^(1/2) + 208*corrincphi*muinc*muphi^4*sig2inc^(1/2)*
        sig2phi^(1/2) + 768*T*corrincphi*muinc^2*muphi^4*sig2inc^(1/2)*
        sig2phi^(1/2) - 1200*corrincphi*muinc*muphi^5*sig2inc^(1/2)*
        sig2phi^(1/2) - 1760*T*corrincphi*muinc^2*muphi^5*sig2inc^(1/2)*
        sig2phi^(1/2) + 4112*corrincphi*muinc*muphi^6*sig2inc^(1/2)*
        sig2phi^(1/2) + 2912*T*corrincphi*muinc^2*muphi^6*sig2inc^(1/2)*
        sig2phi^(1/2) - 9376*corrincphi*muinc*muphi^7*sig2inc^(1/2)*
        sig2phi^(1/2) - 3520*T*corrincphi*muinc^2*muphi^7*sig2inc^(1/2)*
        sig2phi^(1/2) + 14976*corrincphi*muinc*muphi^8*sig2inc^(1/2)*
        sig2phi^(1/2) + 3072*T*corrincphi*muinc^2*muphi^8*sig2inc^(1/2)*
        sig2phi^(1/2) - 17056*corrincphi*muinc*muphi^9*sig2inc^(1/2)*
        sig2phi^(1/2) - 1856*T*corrincphi*muinc^2*muphi^9*sig2inc^(1/2)*
        sig2phi^(1/2) + 13728*corrincphi*muinc*muphi^10*sig2inc^(1/2)*
        sig2phi^(1/2) + 704*T*corrincphi*muinc^2*muphi^10*sig2inc^(1/2)*
        sig2phi^(1/2) - 7488*corrincphi*muinc*muphi^11*sig2inc^(1/2)*
        sig2phi^(1/2) - 128*T*corrincphi*muinc^2*muphi^11*sig2inc^(1/2)*
        sig2phi^(1/2) + 2496*corrincphi*muinc*muphi^12*sig2inc^(1/2)*
        sig2phi^(1/2) - 384*corrincphi*muinc*muphi^13*sig2inc^(1/2)*
        sig2phi^(1/2) + 2*T*muinc^3*sig2phi - 28*T*muinc^3*muphi*sig2phi + 
       12*muinc^2*muphi^2*sig2phi + 172*T*muinc^3*muphi^2*sig2phi - 
       176*muinc^2*muphi^3*sig2phi - 640*T*muinc^3*muphi^3*sig2phi + 
       1084*muinc^2*muphi^4*sig2phi + 1616*T*muinc^3*muphi^4*sig2phi - 
       3832*muinc^2*muphi^5*sig2phi - 2912*T*muinc^3*muphi^5*sig2phi + 
       8888*muinc^2*muphi^6*sig2phi + 3808*T*muinc^3*muphi^6*sig2phi - 
       14464*muinc^2*muphi^7*sig2phi - 3584*T*muinc^3*muphi^7*sig2phi + 
       16936*muinc^2*muphi^8*sig2phi + 2336*T*muinc^3*muphi^8*sig2phi - 
       14160*muinc^2*muphi^9*sig2phi - 960*T*muinc^3*muphi^9*sig2phi + 
       8112*muinc^2*muphi^10*sig2phi + 192*T*muinc^3*muphi^10*sig2phi - 
       2880*muinc^2*muphi^11*sig2phi + 480*muinc^2*muphi^12*sig2phi)/
     (T^2*muinc^4*(1 - 2*muphi + 2*muphi^2)^4))/
  (((-1 + 2/(1 - E^((delta*muinc)/((-1 + muphi)*muphi))))*delta*(1 - muphi)*
       muphi)/T + (((4*E^((2*delta*muinc)/((-1 + muphi)*muphi))*delta^3)/
          ((1 - E^((delta*muinc)/((-1 + muphi)*muphi)))^3*T*(1 - muphi)*
      muphi) + (2*E^((delta*muinc)/((-1 + muphi)*muphi))*delta^3)/
          ((-1 + E^((delta*muinc)/((-1 + muphi)*muphi)))^2*T*(1 - muphi)*
            muphi))*sig2inc)/2 + 
    corrincphi*((-4*E^((2*delta*muinc)/((-1 + muphi)*muphi))*delta^2*
          ((delta*muinc)/((1 - muphi)*muphi^2) - 
            (delta*muinc)/((-1 + muphi)^2*muphi)))/
        ((1 - E^((delta*muinc)/((-1 + muphi)*muphi)))^3*T) - 
       (2*E^((delta*muinc)/((-1 + muphi)*muphi))*delta^2*
          ((delta*muinc)/((1 - muphi)*muphi^2) - 
            (delta*muinc)/((-1 + muphi)^2*muphi)))/
        ((-1 + E^((delta*muinc)/((-1 + muphi)*muphi)))^2*T))*sig2inc^(1/2)*
     sig2phi^(1/2) + sig2phi + (((-2*
            (-1 + 2/(1 - E^((delta*muinc)/((-1 + muphi)*muphi))))*delta)/T + 
         (4*E^((delta*muinc)/((-1 + muphi)*muphi))*delta*
            ((delta*muinc)/((1 - muphi)*muphi^2) - 
              (delta*muinc)/((-1 + muphi)^2*muphi))*(1 - muphi))/
          ((-1 + E^((delta*muinc)/((-1 + muphi)*muphi)))^2*T) - 
         (4*E^((delta*muinc)/((-1 + muphi)*muphi))*delta*
            ((delta*muinc)/((1 - muphi)*muphi^2) - 
              (delta*muinc)/((-1 + muphi)^2*muphi))*muphi)/
          ((-1 + E^((delta*muinc)/((-1 + muphi)*muphi)))^2*T) + 
         (2*E^((delta*muinc)/((-1 + muphi)*muphi))*delta*
            ((-2*delta*muinc)/((1 - muphi)*muphi^3) + 
              (2*delta*muinc)/((-1 + muphi)^2*muphi^2) - 
              (2*delta*muinc)/((1 - muphi)^3*muphi))*(1 - muphi)*muphi)/
          ((-1 + E^((delta*muinc)/((-1 + muphi)*muphi)))^2*T) + 
         (4*E^((2*delta*muinc)/((-1 + muphi)*muphi))*delta*
            ((delta*muinc)/((1 - muphi)*muphi^2) - 
               (delta*muinc)/((-1 + muphi)^2*muphi))^2*(1 - muphi)*muphi)/
          ((1 - E^((delta*muinc)/((-1 + muphi)*muphi)))^3*T) + 
         (2*E^((delta*muinc)/((-1 + muphi)*muphi))*delta*
            ((delta*muinc)/((1 - muphi)*muphi^2) - 
               (delta*muinc)/((-1 + muphi)^2*muphi))^2*(1 - muphi)*muphi)/
          ((-1 + E^((delta*muinc)/((-1 + muphi)*muphi)))^2*T))*sig2phi)/2)

                   

eq5p2Minusinc ::= biaselDivT$pIninc

biaselDivT$pIninc := phi/delta - inc 

eq5p2a ::= biasel01DivT$pIninc

biasel01DivT$pIninc := ((1 - E^((delta*inc)/((-1 + phi)*phi)))*
                       (1 - phi)*phi)/delta - inc 

eq5p4 ::= En$p

En$p := (T*((1 - E^((delta*inc)/(-1 + phi)))*(1 - phi) + phi))/delta

eq5p5 ::= biasnDivsegsInprev

biasnDivsegsInprev := (1 - E^((delta*inc)/(-1 + phi)))*(1 - phi)

eq5p5a ::= biasn11DivsegsInprev

biasn11DivsegsInprev := -((-1 + E^((delta*inc)/(-1 + phi)))^2*
                         (-1 + phi))

eq5p6 ::= biasnDivsegsIninc

biasnDivsegsIninc := -inc + ((1 - E^((delta*inc)/(-1 + 
                      phi)))*(1 - phi) + phi)/delta

eq5p6a ::= biasn01DivsegsIninc

biasn01DivsegsIninc := -inc + ((1 - E^((delta*inc)/(-1 + phi)))*
                       (T/delta - (T*((1 - E^((delta*inc)/(-1 + 
                       phi)))*(1 - phi) + phi))/delta))/T

(*
----------------------------------------------------------------------
      Glossary Adapted for Mathematica Notational Compromises

                      Single Renewal Process

  d0 Duration in On state for qth occurrence (mean mu1)    
  d1 Duration in Off state for qth occurrence (mean mu0)   
  w  = Waiting time between Off-to-On transitions (moments muw,sig2w,mu3w)
  fw  Waiting time probability density function   
  hw  Hazard function corresponding to f(w)       

                      Behavioral Parameters
  inc = 1/muw, Incidence                         
  phi = mu1/(mu1 + mu0), Prevalence                
  mu1  Event Duration                             
  pi   Probability an event is of a specified Type, Pr{"success"}

                      Continuous Observation
  T   Length of observation interval for observing individual p
  N      Count of Off-to-On transitions in interval Tp 
  NDivT   Empirical rate of Off-to-On transitions in interval Tp
  vpT  Variance-time curve, sampling variance of Np(T)/Tp
  L    Time accumulated in On state in interval Tp  
  LDivT  Empirical prevalence in interval Tp        
  Z  Accumulated "successes" in interval Tp        
  ZDivN  Empirical proportion

Multiple Occasions occ; plus subscript indicates a sum over occasions index 
  Npo  Counts in interval Tpo on occasion o         
  s  = {Tp1 , ..., Tpocc},  Observation schedule (1 x occ vector)
  s2  1 x occ vector with elements {(Tp1)2 , ..., (Tpocc)2}
  vpos   A 1 x occ vector with elements vpo(T)          

                           Recorder Errors
  Mpr  Indicates reported count for individual p "Mangled" by errors 
       from recorder r
  scriptR  Denotes a roster of multiple (size R) recorders   
  Mbar Indicates average of reported counts from a roster of recorders
  al   Recorder allocation matrix over multiple occasions (occ x occ)

       Representations for Recorder Errors, Heterogeneity Over Occasions
  kap  Thinning parameter for recorder r error process  
  gam  Superposition parameter for recorder r error process
  om   Imbalance in recorder r error processes
  psi   Proportional variability from recorder r error processes 
  muom, sig2om, mupsi  Moments taken over the population of recorders  
  tau2   variance of  Ipo about Ip.  Distortion of Incidence on occasion o 
         from heterogeneity over occasions  
          
                  Examples of Renewal Processes
  Gamma point process: Index, alpha ; Incidence, inc             
  Alternating Poisson: Incidence, inc ; Prevalence, phi
                       (mu1 = phi/inc , mu0 = (1 - phi)/inc)

     Population of Individual Processes and Psychometric Quantities

  muinc  Mean Incidence over the population of individuals, Ep(Ip)
  sig2inc  Variance of Incidence over population of individuals, Varp(Ip)
  EpvpT    Average sampling variance of Np(T)/Tp over the      
                population of individuals  
  muphi, sig2phi Moments of Prevalence over the population of individuals
  rhoMM'  Inter-recorder correlation; correlation between counts
            reported by recorders r and r'

                    G Theory quantities for pxoxr Design

  porVC = {porVCp, porVCo, porVCr, porVCpo, porVCpr, porVCro,
           porVCpor, porVCeps}                                           
  Variance components for pxoxr design     

  DEL , del  G theory absolute and relative errors

  porVCabs = porVCo + porVCr + porVCpo + porVCpr + porVCro +
             (porVCpor + porVCeps)
  porVCrel = porVCpo + porVCpr + (porVCpor + porVCeps)
  Definitions of absolute and relative variances   


                 Alternatives to Continuous Observation

  delta   Time segment between observations in Instantaneous sampling
        Time segment for recording in One-zero (interval) sampling
  segs  = T/delta, Number of time segments
  el   Count of 1's for spacing delta and observation     
       interval T in Instantaneous sampling 
  el01 Count of 0-to-1 transitions for spacing delta and
       observation interval T in Instantaneous sampling record
  n  Count of 1's in One-zero sampling, recording segment delta and
     observation interval T
  n01 , n11  Counts of 0-to-1 transitions and of 1-to-1 transitions 
             in the One-zero sampling record, recording segment delta 
             and observation interval T

*)
