# SimpleDouble1
# Simple Regression with equivalent double measurement
sem = 'http://www.utstat.toronto.edu/~brunner/openSEM/sage/sem.sage' 
load(sem)
# load('~/sem.sage') For local version: Very handy on the bus.

# Latent regression model
PHIx = ZeroMatrix(1,1); PHIx[0,0]=var('phi'); PHIx
BETA = ZeroMatrix(1,1); BETA[0,0]=var('beta'); BETA 
PSI = ZeroMatrix(1,1); PSI[0,0]=var('psi'); PSI
PHI = RegressionVar(PHIx,BETA,PSI); PHI

# Measurement Model
LAMBDA = ZeroMatrix(4,2)
LAMBDA[0,0] = 1
LAMBDA[1,0] = 1
LAMBDA[2,1] = 1
LAMBDA[3,1] = 1; LAMBDA 
OMEGA = ZeroMatrix(4,4)
OMEGA[0,0] = var('omega11'); OMEGA[1,1] = var('omega11') 
OMEGA[2,2] = var('omega22'); OMEGA[3,3] = var('omega22')
OMEGA[0,2] = var('omega12'); OMEGA[2,0] = var('omega12')
OMEGA[1,3] = var('omega12'); OMEGA[3,1] = var('omega12'); OMEGA
SIGMA = FactorAnalysisVar(LAMBDA,PHI,OMEGA); SIGMA

# Sigma is the covariance matrix under the mis-specified model.
Sigma = SIGMA(omega12=0); Sigma

# detSigma = Sigma.determinant() ; detSigma # Factoring is not helpful
# eigenSigma = Sigma.eigenvalues() # List of eigenvalues
# for item in eigenSigma: show(item)

# Fisher information is the Hessian of the minus log likelihood evaluated at theta.
# So calculate the minus log likelihood, constants and all.
var('n')
detSigma = Sigma.determinant() ; detSigma # Factoring is not helpful
Sigmahat = SymmetricMatrix(4,'s'); Sigmahat
tr = factor((Sigmahat*Sigma.inverse()).trace()); tr
mLogLike = factor( n/2 * log(detSigma) + 4*n/2*log(2*pi) + n/2*tr )
mLogLike

# Differentiate minus log likelihood and set to zero. 
# First need a list of parameters.
param = [phi, beta, psi, omega11, omega22]; param

# Gradient is a vector of partial derivatives
gr = mLogLike.gradient(param) # Gradient
for item in gr: show(factor(item))

# Straight solve never finished.
# Assemble a list of polynomials to give to GroebnerBasis
# Denominators are never zero, and all terms have a factor of n.
poly = [] # Empty list of polynomials
for item in gr:
    niceitem = factor(item)
    term = niceitem * denominator(niceitem)/n
    show(term)
    poly.append(term)
poly[1] = poly[1]/phi  # Second one has a factor of phi

# For Groebner basis,
# Put the s_ij first in the list of variables, unlike what we would do
# with sigma_ij. When I first tried this I did not understand why it was
# such a good choice.
param2 = copy(param) # Work with a copy to avoid changing the original.
sij = Parameters(Sigmahat)
param2.extend(sij) # Put s_ij values at the end
param2.reverse();  # Reversed order of interest
param2

# Try a GroebnerBasis.
basis1 = GroebnerBasis(poly,param2)
len(basis1) # 9

# Take a look
for item in basis1: factor(item)

# Comments
#   Some basis polynomials are huge!
#   The first polynomial has 2 factors. The first (squared) can never be zero in the parameter space. It is possible to see this so easily because s_ij were lowest priority. The second factor says phi=s12.
#   The second  polynomial has 2 factors. The first factor is repeated from polynomial one, and  can never be zero in the parameter space. The second factor says omega11 = (s11+s22)/2 - s12. This is very natural.
#   The third (big) polynomial is a 4th degree polynomial in beta, so in the world of complex numbers it has 4 roots for fixed values of the other variables. Maybe it will simplify if I substitute what we have already found from the first two. 
number3 = Simplify( basis1[2](phi=s12,omega11 = (s11+s22)/2 - s12) )
number3

# Hummm, it remains a 4th degree polynomial, with a few terms involving psi and omega22. Set it aside for now.
# The 4th Groebner basis polynomial is a product of two terms. The first is that same one that appears in the first two polynomials. The second term is promising. Cancel the first term.
print(factor(basis1[0]))

term1 = 2*beta^2*omega11*phi + omega11*omega22 + 2*omega22*phi + 2*omega11*psi + 4*phi*psi
term2 = factor(-basis1[3]/term1) ; term2

term1 = 2*beta^2*omega11*phi + omega11*omega22 + 2*omega22*phi + 2*omega11*psi + 4*phi*psi
term2 = factor(-basis1[3]/term1) ; term2
Simplify( term2(phi=s12,omega11=(s11+s22)/2-s12) )

# That is so beautiful. We now have beta = (s13+s14+s23+s24)/(4*s12), the average.
# The fifth basis polynomial is the product of two terms. The first one is our old friend term1, squared. The second, which thrw me off before because it's directly one of the covariance structure polynomials, says psi = s34 - phi*beta^2
psihat = Simplify( (s34-phi*beta^2)(phi=s12,beta=(s13+s14+s23+s24)/(4*s12)) ); psihat

# So, psihat is not very easy to look at but now it's easy to use. Take the same approach to
# the other MLEs we have so far. 
phihat = s12
omega11hat = (s11+s22)/2 - s12
betahat = (s13+s14+s23+s24)/(4*s12)
# The only parameter for which we do not yet have a unique solution is omega22.
# Try to get it from the third polynomial that was set aside.
Simplify( number3(beta=betahat,psi=psihat) )

# Assuming that variances of observable variables are non-zero, the only solution is
# omega22=0, outside the parameter space. So the third polynomial is out.
# Substitute into the remaining basis polynomials. None of them factors.
# These killing tasks are easy with Sage.
Simplify( basis1[5](phi=phihat,omega11=omega11hat,beta=betahat,psi=psihat) )

# So the sixth one is a quadratic in omega22. The solutions are easy to get but may 
# be unwieldy. Look at the other 3 before trying it.
Simplify( basis1[6](phi=phihat,omega11=omega11hat,beta=betahat,psi=psihat) )

# This will do it. Second term will be non-zero with probability one.
print(_)

omega22hat = (s11*s33 + 2*s12*s33 + s22*s33 - 2*s11*s34 - 4*s12*s34 - 2*s22*s34 + s11*s44 + 2*s12*s44 + s22*s44)/(8*s12)
# Look at that quadratic again (6th)
Simplify( basis1[5](phi=phihat,omega11=omega11hat,beta=betahat,psi=psihat,omega22=omega22hat) )

# Tells me almost nothing. Continue to the end.
Simplify( basis1[7](phi=phihat,omega11=omega11hat,beta=betahat,psi=psihat) )

# Ah, this one is so much simpler and nicer, corresponding to the smart way of identifying omega22. The first omega22hat might not even be in the parameter space. I could let n -> infinity to see.
target = SigmaOfTheta(Sigma,'s'); target
Simplify(omega22hat(target))

# Typo? Seventh should be zero at these MLEs.
Simplify( basis1[6](phi=phihat,omega11=omega11hat,beta=betahat,psi=psihat,omega22=omega22hat) )

# Well, it's zero. How educational. Since MLEs are consistent, I guess it means that the supposed omega22hat was not an MLE. Maybe it was a minimum or saddle point. Could check this later. Do I really need to verify that my putative MLEs are consistent under the model? Better do it.
omega11hat(target)
betahat(target)
# Those were obvious. Psihat actually requires a check.
Simplify( psihat(target) ) # It had me scared for a moment without Simplify.

# Take a better omega22hat, but save the old one as a specimen.
badomega22hat = copy(omega22hat); badomega22hat
omega22hat = (s33+s44)/2 - s34
Simplify( basis1[7](phi=phihat,omega11=omega11hat,beta=betahat,psi=psihat,omega22=omega22hat) )

# Last one is zero
Simplify( basis1[8](phi=phihat,omega11=omega11hat,beta=betahat,psi=psihat,omega22=omega22hat) )

# Go back to the quadratic 6th again, with the right omega22hat.
Simplify( basis1[5](phi=phihat,omega11=omega11hat,beta=betahat,psi=psihat,omega22=omega22hat) )

# Non-zero with probability one -- another phantom solution?
# Verify that the MLEs really do make the derivatives exactly zero.
for item in gr: Simplify( item(phi=phihat,omega11=omega11hat,beta=betahat,psi=psihat,omega22=omega22hat) )

# Have consistent estimators that make the derivatives of the log likelihood function all equal zero.
# These must be the MLEs. List them nicely so I don't need to re-compute this.
# Thetahat will be in same order as param
thetahat = [phihat, betahat, psihat, omega11hat, omega22hat]
for j in range(5): param[j],thetahat[j]

for j in range(5): 
    print(latex(thetahat[j]))
    print(' ')

# Where are the MLEs going under the true model?
truetarget = SigmaOfTheta(SIGMA,'s')
for j in range(5): param[j],Simplify(thetahat[j](truetarget))

# Calculate the reproduced covariance matrix
SigmaOfThetaHat = factor( Sigma(phi=phihat, beta=betahat, psi=psihat, omega11=omega11hat, omega22=omega22hat) )
SigmaOfThetaHat

# Large-sample target of SigmaOfThetaHat under the true model
SigmaOfThetaHatTarget = SigmaOfThetaHat(truetarget)
SigmaOfThetaHatTarget

# Where is fmax going? Everything will depend on the product of SigmaHat times the inverse of SigmaOfThetaHat. It's likely to be ugly, but we can hope for simplicity in the limit.
Prod = factor( SIGMA*SigmaOfThetaHatTarget.inverse() )
Prod

# I believe I am having a religious experience.
Prod.determinant()

# So the large-sample target of fmax is -log(1-corr^2), where corr is the correlation between e1 and e2 (also e3 and e4). In the SAS output, had RMSEA = 0.2275 with corr=0.50
sqrt(-log(1-0.5^2)/5 - 1/999) # 0.239867 

# Cohen's (1988) definition of a "small" correlation is r = 0.1
sqrt(-log(1-0.1^2)/5) # 0.0448

# So for this problem, RMSEA < 0.044 could be an indication of adequate fit.

# What is the power of the goodness of fit test corr = 0.1 and n = 500?
# Noncentrality parameter is n times the target of fmax, which is 
# -log(1-corr^2).
corr = omega12/sqrt(omega11*omega22)
fmaxinf = -log(1-corr^2); fmaxinf
ncp = n*fmaxinf
ncp

nncp = ncp(n=500, omega11=1, omega22=3, omega12 = 0.1732050808)
nncp # 5.02516792926493

# SageMath does not seem to have the non-central chi-squared distribution.
# Invoke R from within SageMath. Using %r runs a little R session in that cell.

%r 
corr = 0.1; n = 500
fmaxinf = -log(1-corr^2); ncp = n*fmaxinf; ncp
critval = qchisq(0.95,5)
1 - pchisq(critval,5,ncp) # Asymptotic power of the goodness of fit test = 0.3644613
 
# The r() command produces output in an "R data structure" that cannot be viewed when the Typeset option is checked. The error message says that the R package Hmisc is needed, but my attempts to install it in the SageMath environment were unsuccessful, and painful. However, you can convert an R-type object back to Python/Sage. An example follows. 
r('rnorm(5)')

print(r('rnorm(5)')) # Plain text output

show(corr) # Symbolic Sage version
print(r('corr')) # R numerical version, defined earlier

# So apparently the R session is still running in the background! Get the Asymptotic 
# power into the SageMath environment.
apow = r('1 - pchisq(critval,5,ncp)') # An R data structure
print(apow)

Apower = apow._sage_(); Apower

# That was excellent. Now how can one get SageMath numbers into R? Recall that 
# nncp is a numerical version of ncp. Convert it to a string with str(nncp). 
# The r() command accepts string as input. 
r('NCP = '+str(nncp)) # An R assignment statement
print(r('NCP'))

print( r('1 - pchisq(critval,5,NCP)') ) # Using NCP, which is now defined in the R session.

# Numerical version of the true covariance matrix, for input to SAS. Let corr=0.1.
# omega12 = corr*sqrt(omega11*omega22) = 0.1*sqrt(1*3) = 0.1732051
nSIGMA = SIGMA(beta=0,phi=2,omega11=1,omega22=3,omega12=0.1732051,psi=1)
nSIGMA

print(nSIGMA) # To edit into SAS input

print(latex(nSIGMA))

# Numerical versions of MLEs
for j in range(5): param[j],Simplify(thetahat[j](truetarget)), thetahat[j](truetarget)(beta=0,phi=2,omega11=1,omega22=3,omega12=0.1732051,psi=1)

# Strictly numerical version of the large-sample targets
for j in range(5): 
    print(thetahat[j](truetarget)(beta=0,phi=2,omega11=1,omega22=3,omega12=0.1732051,psi=1))