######################################################################
##SameSexMarriages: Save this file as  SameSexMarriages              #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read SameSexMarriages <Enter>                      #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University,                   #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: April 27, 2011
 
print(`Created:  April 27, 2011`):
print(` This is SameSexMarriages, `):
print(`it is a Maple program that generated the article `):
print(`The Number of Same-Sex Marriages in a Perfectly Bisexual`):
print(`Population is Asymptotically Normal`):
print(`by Shalosh B. Ekhad`):
print(`published in the prestigious on-line journal`):
print(`the Personal Journal of Shalosh B. Ekhad and Doron Zeilberger`):
print(``):
print(`The URL of this Maple package  is:`):
print(`http://www.math.rutgers.edu/~zeilberg/tokhniot/SameSexMarriages`):
print(``):
print(`It requires the packages HISTABRUT and GuessRat:`):
print(`http://www.math.rutgers.edu/~zeilberg/tokhniot/HISTABRUT`):
print(`http://www.math.rutgers.edu/~zeilberg/tokhniot/GuessRat`):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):
 print(`For a list of the procedures type Help();, for help with`):
 print(`a specific procedure, type Help(procedure_name);   .`):
 print(``):

with(combinat):

read HISTABRUT:
read GuessRat:

Help:=proc()

if args=NULL then
 print(`The main procedures are: `):
 print(`F, SSM, SSMlatex `):


elif nops([args])=1 and op(1,[args])=F then

print(`F(n,x): the probability generating function for the random`):
print(`variable "number of same-sex marriages" in  a perfectly`):
print(`bi-sexual population of 2n men and 2n women`):
print(`For example, try:`):
print(`F(10,x);`):

elif nops([args])=1 and op(1,[args])=SSM then
print(`SSM(n,L,M,K): Inputs a symbol n and a positive `):
print(`integers L,M, and K`):
print(`that generates a mathematics article that`):
print(`semi-rigrously proved the asymptotic normality`):
print(`of the random variable "number of same-sex marriages`):
print(`in a perfectly bi-sexual population, by going up`):
print(`to the L-th moment, and displaying asymptotic expansions`):
print(`of these normalized moments up to order M.`):
print(`K is the size of the data used for guessing.`):
print(`It also outputs these explicit expressions`):
print(`for the mean, variance, and normalized moments`):
print(`For example, try:`):
print(`SSM(n,2,4,30);`):


elif nops([args])=1 and op(1,[args])=SSMlatex then
print(`SSMlatex(n,L,M,K): like SSM(n,L,M,K) but prints it out`):
print(`in latex, so it could be copied-and-pasted into the real article`):
print(`For example, try:`):
print(`SSMlatex(n,2,4,30);`):


else
print(`There is no ezra for`,args):
fi:
 
end:


#F(n,x): the probability generating function for the random
#variable "number of same-sex marriages" in  a perfectly
#bi-sexual population of 2n men and 2n women
#For example, try:
#F(10,x);
F:=proc(n,x) local k,gu,lu:
gu:=add(binomial(2*n,2*k)^2*(2*n-2*k)!*((2*k)!/k!/2^k)^2*x^(2*k),k=0..n):
lu:=(4*n)!/(2*n)!/2^(2*n):
gu/lu:
end:

#SSM(n,L,M,K): Inputs a symbol n and a positive 
#integers L,M, and K
#that generates a mathematics article that
#semi-rigrously proved the asymptotic normality
#of the random variable "number of same-sex marriages"
#in a perfectly bi-sexual population, by going up
#to the L-th moment, and displaying asymptotic expansions
#of these normalized moments up to order M.
#K is the size of the data used for guessing.
#It also outputs these explicit expressions
#for the mean, variance, and normalized moments
SSM:=proc(n,L,M,K) local gu,mu,lu,i,x,mu1,mu2,sof,j,as1,k,t0:
t0:=time():
mu:=[seq(F(i,x),i=1..K)]:

gu:=[seq(AlphaSeq(mu[i],x,L),i=1..K)]:

lu:=[seq(gu[i][1],i=1..K)]:
mu1:=GuessRat(i->lu[i],[n],10,1);

if mu1=FAIL then
print(`Not enough data, make the last argument, K, larger`):
 RETURN(FAIL):
fi:
lu:=[seq(gu[i][2],i=1..K)]:
mu2:=GuessRat(i->lu[i],[n],10,1);
if mu2=FAIL then
print(`Not enough data, make the last argument, K, larger`):
 RETURN(FAIL):
fi:
sof:=[mu1,mu2]:
print(`The Number of Same-Sex Marriages in a Perfectly Bisexual`):
print(`Population is Asymptotically Normal`):
print(``):
print(`by Shalosh B. Ekhad`):
print(``):
print(`Theorem: Consider  a population of 2n men and 2n women`):
print(`where every individual is equally attracted to either sex`):
print(`and chooses his or her mate according to other criteria.`):
print(`Also assume that everyone gets married.`):
print(`Then the expectation of the random variable `):
print(` "Number of same-sex marriages" is`):
print(mu1):
print(`That asymptotically is`):
print(asympt(mu1,n,M)):
print(`And its variance is`):
print(mu2):
print(`That asymptotically is`):
print(asympt(mu2,n,M)):
print(`Furthermore, this random variable is asymptotically normal`):
print(`at least  by looking up to the`, L, `-th normalized moments.`):
print(``):
print(`Semi-Rigorous Proof:`):
print(`The probability generating function is (why?)`):
print(Sum(((4*n)!/(2*n)!/2^(2*n))^(-1)*
binomial(2*n,2*k)^2*(2*n-2*k)!*((2*k)!/k!/2^k)^2*x^(2*k),k=0..n)):
print(`By taking successive derivatives, into the special cases`):
print(`n from`, 1 , `to `, K):
print(`and fitting the data by rational functions, it emerges that the`):
print(`mean and variance are indeed as stated by the theorem`):
print(`as for the higher moments, we have the following`):

for j from 3 to L do 
lu:=[seq(gu[i][j],i=1..K)]:

if j mod 2=1 then
lu:=[seq(lu[i]^2,i=1..K)]:
mu:=GuessRat(i->lu[i],[n],K/2,1);

if mu=FAIL then
print(`Not enough data, make the last argument, K, larger`):
 RETURN(FAIL):
fi:

print(`The normalized`, j, `-th moment about the mean is`):
print(sqrt(mu)):

as1:=asympt(sqrt(mu),n,M):

print(`The asymptotics of this is`):
print(as1):
print(`Notice that it goes to 0, as it should since the odd moments of`):
print(`the normal distribution are all 0`):

sof:=[op(sof),sqrt(mu)]:
elif j mod 2=0 then
mu:=GuessRat(i->lu[i],[n],K/2,1);

if mu=FAIL then
print(`Not enough data, make the last argument, K, larger`):
 RETURN(FAIL):
fi:

print(`The normalized`, j, `-th moment about the mean is`):
print(mu):

as1:=asympt(mu,n,M):
print(`The asymptotics of this is`):
print(as1):

print(`Notice that it converges to`, (j)!/(j/2)!/2^(j/2), `as it should`):
print(` this being the `, j, `-th moment of the normal distibution`):
sof:=[op(sof),mu]:
else
 RETURN(FAIL):
fi:
od:

print(`This took`, time()-t0, `seconds . `):
sof:
end:


#SSMlatex(n,L,M,K): LaTeX version
#Inputs a symbol n and a positive 
#integers L,M, and K
#that generates a mathematics article that
#semi-rigrously proved the asymptotic normality
#of the random variable "number of same-sex marriages"
#in a perfectly bi-sexual population, by going up
#to the L-th moment, and displaying asymptotic expansions
#of these normalized moments up to order M.
#K is the size of the data used for guessing.
#It also outputs these explicit expressions
#for the mean, variance, and normalized moments
SSMlatex:=proc(n,L,M,K) local gu,mu,lu,i,x,mu1,mu2,sof,j,as1,k,t0:
t0:=time():
mu:=[seq(F(i,x),i=1..K)]:

gu:=[seq(AlphaSeq(mu[i],x,L),i=1..K)]:

lu:=[seq(gu[i][1],i=1..K)]:
mu1:=GuessRat(i->lu[i],[n],10,1);

if mu1=FAIL then
 RETURN(FAIL):
fi:
lu:=[seq(gu[i][2],i=1..K)]:
mu2:=GuessRat(i->lu[i],[n],10,1);
if mu2=FAIL then
 RETURN(FAIL):
fi:
sof:=[mu1,mu2]:
print(`The Number of Same-Sex Marriages in a Perfectly Bisexual`):
print(`Population is Asymptotically Normal`):
print(``):
print(`by Shalosh B. Ekhad`):
print(``):
print(`Theorem: Consider  a population of 2n men and 2n women`):
print(`where every individual is equally attracted to either sex`):
print(`and chooses his or her mate according to other criteria.`):
print(`Also assume that everyone gets married.`):
print(`Then the expectation of the random variable `):
print(` "Number of same-sex marriages" is`):
print(latex(mu1)):
print(`That asymptotically is`):
print(latex(asympt(mu1,n,M))):
print(`And its variance is`):
print(latex(mu2)):
print(`That asymptotically is`):
print(latex(asympt(mu2,n,M))):
print(`Furthermore, this random variable is asymptotically normal`):
print(`at least  by looking up to the`, L, `-th normalized moments.`):
print(``):
print(`Semi-Rigorous Proof:`):
print(`The probability generating function is (why?)`):
print(latex(Sum(((4*n)!/(2*n)!/2^(2*n))^(-1)*
binomial(2*n,2*k)^2*(2*n-2*k)!*((2*k)!/k!/2^k)^2*x^(2*k),k=0..n))):
print(`By taking successive derivatives, into the special cases`):
print(`n from`, 1 , `to `, K):
print(`and fitting the data by rational functions, it emerges that the`):
print(`mean and variance are indeed as stated by the theorem`):
print(`as for the higher moments, we have the following`):

for j from 3 to L do 
lu:=[seq(gu[i][j],i=1..K)]:

if j mod 2=1 then
lu:=[seq(lu[i]^2,i=1..K)]:
mu:=GuessRat(i->lu[i],[n],K/2,1);

if mu=FAIL then
 RETURN(FAIL):
fi:

print(`The normalized`, j, `-th moment about the mean is`):
print(latex(sqrt(mu))):

as1:=asympt(sqrt(mu),n,M):

print(`The asymptotics of this is`):
print(latex(as1)):
print(`Notice that it goes to 0, as it should since the odd moments of`):
print(`the normal distribution are all 0`):

sof:=[op(sof),mu]:
elif j mod 2=0 then
mu:=GuessRat(i->lu[i],[n],K/2,1);

if mu=FAIL then
 RETURN(FAIL):
fi:

print(`The normalized`, j, `-th moment about the mean is`):
print(latex(mu)):

as1:=asympt(mu,n,M):
print(`The asymptotics of this is`):
print(latex(as1)):

print(`Notice that it converges to`, latex((j)!/(j/2)!/2^(j/2)), `as it should`):
print(` this being the `, j, `-th moment of the normal distibution`):
sof:=[op(sof),mu]:
else
 RETURN(FAIL):
fi:
od:

print(`This took`, time()-t0, `seconds . `):
sof:
end: