################ A little sage demonstration ################

notebook()
# Select new worksheet option, check typeset

1+1
1.0 + 1.0
1 + 1.0

# Of 100 graduating students, how many ways are there for 60 to be employed in a job related to their field of study, 30 to be employed in a job unrelated to their field of study, and 10 unemployed?

factorial(100)/(factorial(60)*factorial(30)*factorial(10))


# In R,
# prod(1:100)/(prod(1:60)*prod(1:30)*prod(1:10))
# > prod(1:100)/(prod(1:60)*prod(1:30)*prod(1:10))
# [1] 1.165214e+37


f(x) = exp(-x)
integrate(f(x),x)
integrate(f(x),x,0,infinity)

g(x) = 1/(1+x^2) ; integrate(g(x),x)

h(x,t) = 1/t * exp(-x/t) ; integrate(h(x,t),x)

integrate(f(y),y,0,infinity)
var('y')
integrate(f(y),y,0,infinity)

g(x) = exp(x*t) * f(x)
assume(t-1<0)
mgf(t) = integrate(g(y),y,0,infinity) ; mgf(t)

mgf(t) = integrate(g(y),y,0,infinity) 

derivative(mgf(t),t)

derivative(mgf(t),t)(t=0)
derivative(mgf(t),t,2)(t=0) # Second derivative

L = var('lambda') 
# lambda has a special meaning. But if you assign 
# it to a variable you can use it as a symbol.
p(x) = exp(-L) * L^x / factorial(x)
p(x).sum(x,0,oo)
(x*p(x)).sum(x,0,oo) # Expected value


# Matrices...

mat = matrix(SR,2,2); mat # SR means Symbolic Ring

var('x y')
mat[0,0]=x; mat[0,1]=x
mat[1,0] = y; mat[1,1]=y ; mat

# Just paste some in to test

mat.is_square() # T-F test
mat.is_symmetric() # T-F test
mat.det()
mat.inverse() # Error for this ex
mat.is_invertible() # T-F test
mat.characteristic_polynomial()
mat.eigenvalues()
mat.rank()
mat.rows()

mat.diagonal() # No such attributite
# Get the main diagonal
diag = []
for i in srange(0,mat.nrows()): diag.append(mat[i,i])

# Finding max-min
var('x1 x2 mu1 mu2')
f(x1,x2) = (x1-mu1)^2 + (x2-mu2)^2
d1 = derivative(f(x1,x2),x1); d1
d2 = derivative(f(x1,x2),x2); d2
eqns = [d1==0, d2==0] # A list of equations to solve
solve(eqns,[x1,x2]) # Solve for x1 and x2
H = matrix(SR,2,2)
H[0,0] = derivative(d1,x1); H[0,1] = derivative(d1,x2)
H[1,0] = derivative(d2,x1); H[1,1] = derivative(d2,x2)
H; H.determinant() # Positive, so concave up (Rule is necessary but not sufficient if more than 2x2)
# Or an easier way ...
H2 = f.hessian(); H2
H2.determinant()

# Variance of geometric ( a nasty chore)
var('t theta')
mgf(t) = theta*exp(t)/(1-exp(t)*(1-theta)) # By hand, unfortunately
m1 = derivative(mgf,t); m1
m1 = factor(expand(m1)); m1
m2 = derivative(mgf,t,2); m2
m2 = factor(expand(m2)); m2
m2(t=0)
v = m2(t=0)-(m1(t=0))^2; v
factor(expand(v))

# Developing the function ml
# Given a likelihood function and a parameter vector, return the MLE as a dictionary.

def MLE(like,params):
    """
    Given a likelihood function, return the presumed MLE in the form of a 
    dictionary. I say presumed because all it does is diffferentiate the minus 
    log likelihood function, set the derivatives to zero and solve. There is no 
    attempt to check whether it's a minimum. Later versions of the function may 
    allow this.
    
    Symbolic variables must be declared, and the function like must give the 
    parameters explicitly as arguments, as in the example below. Notice the 
    symbols for sufficient statistics. As far as I can tell this step must be 
    done by hand; apparently, Sage can't do it.
    
    Example:
    
    var('theta1 theta2 x1 x2 n')
    # Multinomial
    mlike(theta1,theta2) = theta1^x1 * theta2^x2 * (1-theta1-theta2)^(n-x1-x2)
    theta = [theta1, theta2]
    multimle = MLE(mlike,theta); multimle

    """
    ell = -log(like)
    grad = ell.gradient()
    eqns = [] # Empty list
    for a in grad: eqns.append(a==0)
    MLE = solve(eqns,params,solution_dict=True)[0]
    # Solve returns a list of dictionaries. This one (the only one) 
    # is number zero.
    return MLE


var('theta1 theta2 x1 x2 n')
mlike(theta1,theta2) = theta1^x1 * theta2^x2 * (1-theta1-theta2)^(n-x1-x2)
theta = [theta1, theta2]
MLE(mlike,theta)
loglike.simplify_full()


H = matrix(SR,2,list(H)) # Symbolic ring version of H. Can evaluate at a dictionary.
Hmle = factor(expand(H(mle))); Hmle


def ml(like,params):
    """
    Returns the mle in the form of a dictionary, the Hessian, and the inverse
    of the Hessian evaluated at the MLE, which is the estimated asymptotic 
    covariance matrix of the MLE.
    """
    ell = -log(like)
    grad = ell.gradient()
    eqns = [] # Empty list
    for a in grad: eqns.append(a==0)
    mle = solve(eqns,params,solution_dict=True)[0]
    # Solve returns a list of dictionaries. This one (the only one) 
    # is number zero.
    k = len(params)
    H = ell.hessian()
    H = matrix(SR,k,list(H)) # Symbolic ring version of H. Can 
                             # evaluate at a dictionary.
    H = factor(expand(H))
    Hmle = H(mle)
    Vhat = factor(expand(Hmle.inverse()))
    ml = [mle,H,Vhat]
    return ml

var('theta1 theta2 x1 x2 n')
mlike(theta1,theta2) = theta1^x1 * theta2^x2 * (1-theta1-theta2)^(n-x1-x2)
theta = [theta1, theta2]
max = ml(mlike,theta); max
Hess = max[1]; Hess
I = Hess(x1 = n*theta1,x2=n*theta2)
I = factor(expand(I)); I
V = factor(expand(I.inverse()))

# Normal
var('tau,sigma,mu,m,s,pi,n')
nlike(tau,mu) = tau^(-n/2) * (2*pi)^(-n/2) * exp(-n/2 * (s+(m-mu)^2)/tau)
theta = [mu,tau]
MLE(nlike,theta)
likeit = ml(nlike,theta); likeit
H = likeit[1]
H = H(tau = sigma^2); H
H(m^2=sigma^2+mu^2) # Nope
print(H[0,0])
print(H[1,1])
L = [ [1/2*(2*(sigma^2+mu^2) - 4*mu*mu + 2*mu^2 - sigma^2 + 2*sigma^2)*n/sigma^6,0],
[0,n/sigma^2] ] # List of lists
I = matrix(SR,L); I