# HW 20
# Jike Liu


read C20.txt;


# Question 1: 

# GenBSsamp(L,k): a bootstrap sample of length k*nops(L)
# Assumes that k*nops(L) is a positive integer
GenBSsamp:=proc(L,k) local n,m,ra,i:
n:=nops(L):
m:=k*n:

if not type(m,posint) then
   error "k*nops(L) must be a positive integer":
fi:

ra:=rand(1..n):
[seq(L[ra()],i=1..m)]:
end:


# GenBSse(L,B,k): bootstrap estimate of the standard error
# using B bootstrap samples, each of length k*nops(L)
GenBSse:=proc(L,B,k) local M,mu,i:
M:=[seq(Av(GenBSsamp(L,k)),i=1..B)]:
mu:=Av(M):
evalf(sqrt(add((M[i]-mu)^2,i=1..B)/(B-1))):
end:

# Test:
L:=CH():
SE(L);
# 2.635231383

GenBSse(L, 1000, 1/2);
# 3.733414128

GenBSse(L, 1000, 1);
# 2.657853388

GenBSse(L, 1000, 2);
# 1.869211489

# GenBSse(L,1000,1) is closest to SE(L)


# Question 2:

with(combinat):

# maj(pi): the major index
maj:=proc(pi) local n,i,co:
co:=0:
n:=nops(pi):
for i from 1 to n-1 do
 if pi[i]>pi[i+1] then
   co:=co+i:
 fi:
od:
co:
end:

# inv(pi): number of inversions
inv:=proc(pi) local n,i,j,co:
n:=nops(pi):
co:=0:
for i from 1 to n do
 for j from i+1 to n do
   if pi[i]>pi[j] then
      co:=co+1:
   fi:
 od:
od:
co:
end:

# SampInvMaj(n,K): sample of inv(pi)*maj(pi) over K random permutations of length n
SampInvMaj:=proc(n,K) local i,pi,L:
L:=[]:
for i from 1 to K do
   pi:=randperm(n):
   L:=[op(L), inv(pi)*maj(pi)]:
od:
L:
end:

L:=SampInvMaj(100,1000):
# average
Av(L);
# standard error
SE(L);
# 6.166902469*10^6
# 19682.76683


# E(inv(pi)*maj(pi)) = (binomial(n,2)*(binomial(n,2)+1))/4.

# Proof:
# Let N = binomial(n,2) = n*(n-1)/2.

# We use the identity

# E(inv(pi)*maj(pi)) = Cov(inv(pi),maj(pi)) + E(inv(pi))*E(maj(pi)).

# In: https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimPDF/invmaj3rd.pdf
# Baxter-Zeilberger gives the exact covariance: 

# Cov(inv(pi),maj(pi)) = n*(n-1)/8 = N/4.

# Also, inv and maj are equidistributed on S_n, so they have the same mean. Since

# E(inv(pi)) = n*(n-1)/4 = N/2,

# it follows that

# E(maj(pi)) = n*(n-1)/4 = N/2.

# Therefore

# E(inv(pi)*maj(pi))
# = Cov(inv(pi),maj(pi)) + E(inv(pi))*E(maj(pi))
# = N/4 + (N/2)*(N/2)
# = N/4 + N^2/4
# = N*(N+1)/4.

# Substituting N = binomial(n,2), we get

# E(inv(pi)*maj(pi)) = (binomial(n,2)*(binomial(n,2)+1))/4.

# This completes the proof.