######################################################################
##GANAVIM: Save this file as  GANAVIM                                #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read GANAVIM<Enter>                                #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: alef-bet tishrei 5773, while listening to the prayers in
#print(`http://www.math.rutgers.edu/~zeilberg/noraim.html  .`):
 
print(`Created:  alef-bet tishrei 5773 (heh-tav-shin-ayin-gimel)`):
print(` This is GANAVIM `):
print(`a Maple package that uses symbol-crunching`):
print(`for Statistical Analysis of Empirical Date`):
print(`(Illustrated by Analyzing the immoral practice of the Hillel Foundation to Charge Admission to High Holiday Services from Non-Students)`):
print(``):
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 ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are:  Etnan, GenData, NoHillel`):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: GF1, StatAnal, StatAnalA `):
 print(`SipurStatAnal, SipurStatAnalA `):
 print(` `):


elif nops([args])=1 and op(1,[args])=Etnan then
print(`Etnan(NA): The data-list of Hillel's that have services open to the public (either free or for a charge)`):

elif nops([args])=1 and op(1,[args])=GenData then
print(`GenData(): the empirically derived data of the "top" fifty colleges according`):
print(`to the notorious US News and World ranking. It is a list of length 50 each entry consisting of`):
print(`of a list of length 5, [Ranking,Name,Tuition,NumberOfStudents,AmountChargedByHillelToPrayFromNonStudents]`):
print(`try:`):
print(` GenData(); `):

elif nops([args])=1 and op(1,[args])=GF1 then
print(`GF1(L,x): given a list of data, of the same length, all the format [Name,Data1,..., DataN] all of`):
print(`the same length where Data1, ..., DataN are numeric, finds the generating function`):
print(`as a function of the indexed variables x[1], ..., x[N]`):
print(`Try: GF1(Etnan(NA),x);`):


elif nops([args])=1 and op(1,[args])=NoHillel then
print(`NoHillel(L): given a list of length-5-lists returns the list of length-4 lists obtained by`):
print(`deleting the last component. Try: `):
print(`NoHilel(GenData());`):


elif nops([args])=1 and op(1,[args])=SipurStatAnal then
print(`SipurStatAnal(L,r,A,Koteret,POPULATION):`):
print(`verbose version of StatAnal(L,r)`):
print(`where A is a list of attributes (in English)`):
print(`and r is a positive integer.`):
print(`Try `):
print(`SipurStatAnal(NoHillel(GenData()),10,[Ranking,Tuition,NumberOfStudet],StatisticalAnalysisOfTuitionCosts,UScolleges):`):

elif nops([args])=1 and op(1,[args])=SipurStatAnalA then
print(`SipurStatAnalA(L,r,A,Koteret,POPULATION):`):
print(`verbose version of StatAnalA(L,r)`):
print(`where A is a list of attributes (in English)`):
print(`and r is a positive integer.`):
print(`Try `):
print(`SipurStatAnalA(Etnan(NA),10,[Ranking,Tuition,NumberOfStudents,TicketCostForPraying],StatisticalAnalysisOfTheHillelFoundationTicketCost,UScolleges):`):

elif nops([args])=1 and op(1,[args])=StatAnal then
print(`StatAnal(L,r):Inputs a list of lists L, all of the same length`):
print(`where each entry starts with the name of the individual and`):
print(`numerical values of the same attibutes`):
print(`Outputs:(1) A list of length nops(L)-1 each of whose entries`):
print(`contains [ave,s.d.,alpha_3, ..., alpha_r]`):
print(`for the i-th attribute`):
print(`(2) A list of lists, let's call it M, such that M[i][j]`):
print(`is the correlation between the i-th and j-th attribute`):
print(`Try:`):
print(`StatAnal(NoHillel(GenData()),10);`):

elif nops([args])=1 and op(1,[args])=StatAnalA then
print(`StatAnalA(L,r):Inputs a list of lists L, all of the same length`):
print(`where each entry starts with the name of the individual and`):
print(`numerical values of the same attibutes, where the last one is`):
print(`the most important.`):
print(`Outputs:(1) A list of length r [ave,s.d.,alpha_3, ..., alpha_r]`):
print(`(2) A list of length nops(L[1])-2 consisting of the correlation`):
print(`of the last, important attribute with all the rest`):
print(`Try:`):
print(`StatAnalA(Etnan(NA),10);`):


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



GenData:=proc()
[
[1,Harvard,39849,19627,300],
[1,Princeton,37000,7802,180],
[3,Yale,40500,11701,154],
[4,Columbia,45290,22283,100],
[5,CalTech,37704,2175,NA],
[5,MIT,40732,10566,NA],
[5,Stanford,40569,19535,90],
[5,Chicago,42783,12781,150],
[5,Penn,42098,19842,180],
[10,Duke,41958,14983,100], #suggested
[11,Dartmouth,42996,6141,NA],
[12,NorthWestern,41983,19389,250],
[13,JohnHopkins,42280,21092,NA],
[14,WashingtonUstLewis,41992,13820,0], #appreciated
[15,Brown,42230,8695,NA], #secret
[15,Cornell,41541,20939,NA],
[17,Rice,35551,5879,NA], #Houston
[17,Vanderbilt,41332,12714,NA],
[19,NotreDame,41417,11992,NA],
[20,Emory,41164,13381,150],
[21,UCB,34645,35838,NA], #in-state tuition=11767
[22,Georgetown,41398,16871,NA], #students pay 20
[23,CarnegieMellon,43396,11618,NA],
[23,USC,42818,36896,NA], #$18/service for students!
[25,UCLA,34482,39593,300], #in-stae tuition=11604
[25,Virginia,36570,24391,0], #ask for donations #in-state=11576
[25,WakeForest,41576,7162,NA],
[28,Michigan,37265,41924,NA], #in-state 12590
[29, Tufts,42962,10489,NA], #only for students
[29,UNC,26834,29390,90], #suggested
[31,BostonCollege,42204,14015,NA],
[31, Brandeis,41860,5642,300], #min! suggested=600
[33,WilliamMary,35409,8000,NA], #in-state=13132
[33,NYU,41606,43797,85*2], #95*2 at the door (pay separately)
[35, Rochester,41826,10111,NA],
[36, GeorgiaTech,27862,20720,150], #in Emory, in-state=9652
[37,UCSD, 35006, 29176, 2*75], #separate, in-state=12128
[38, CaseWestern, 39120, 9837,NA],
[38, Lehigh, 40960, 7051, NA],
[38,UCD,35672,31392,NA], #in-state=12794
[38, Miami, 39654,15657,NA],
[42,UCSB,35386,22218,NA], #in-state=12508
[43,UWASeattle,28054, 42446, 100], #in-state=10574 #suggested
[44, Wisconsin, 25421, 42596,NA], #in-state=9671
[45,PennState,28066,45233,0],#in-state=15984, #optional
[45,UCI, 35780,26994,NA], #in-state=12902, #Orange-county Hillel
[45, UIUC, 27700, 43862,0], #in-state=13558
[45, UTexas, 32506, 51195, 180], #in-state=9794
[45, Yeshiva,35150, 6445, NA],
[50,GWU, 42425,25135,125]
]:
end:

#Etnan(NA): The data-list of Hillel's that have services open to the public (either free for a charge)
Etnan:=proc(NA) local gu,mu,i:
gu:=GenData():
mu:=[]:

for i from 1 to nops(gu) do
 if gu[i][5]<>NA then
  mu:=[op(mu),[gu[i][2],gu[i][1],gu[i][3],gu[i][4],gu[i][5]]]:
 fi:
od:

mu:

end:





#NoHillel(L): given a list of length-5-lists returns the list of length-4 lists obtained by
#deleting the last component. Try: NoHilel(GenData());

NoHillel:=proc(L) local i:
[seq([L[i][2],L[i][1],L[i][3],L[i][4]],i=1..nops(L))]:
end:


#GF1(L,x): given a list of data, of the same length, all the format [Name,Data1,..., DataN] all of
#the same length where Data1, ..., DataN are numeric, finds the generating function
#as a function of the indexed variables x[1], ..., x[N]
GF1:=proc(L,x) local i,j:
if nops({seq(nops(L[i]),i=1..nops(L))})<>1 then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if {seq(seq(type(L[i][j],numeric),j=2..nops(L[i])),i=1..nops(L))}<>{true} then
 print(`Bad input`):
 RETURN(FAIL):
fi:

add(mul(x[j-1]^L[i][j],j=2..nops(L[i])),i=1..nops(L)):

end:



#StatAnalA(L,r):Inputs a list of lists L, all of the same length
#where each entry starts with the name of the individual and
#numerical values of the same attibutes, where the last one is
#the most important.
#Outputs:(1) A list of length r [ave,s.d.,alpha_3, ..., alpha_r]
#(2) A list of length nops(L[1])-2 consisting of the correlation
#of the last, important attribute with all the rest
StatAnalA:=proc(L,r) local gu,x,N,zonot,ave,lu,i,sd,r1,Ave,SD,
gu1,menuvalim,kak,i1:

if nops({seq(nops(L[i]),i=1..nops(L))})<>1 then
 print(`Bad input`):
 RETURN(FAIL):
fi:

gu:=GF1(L,x):

N:=nops(L[1])-1:

lu:=subs({seq(x[i]=1,i=1..N-1)},gu):

if subs(x[N]=1,lu)=0 then
 RETURN(FAIL):
fi:

lu:=lu/subs(x[N]=1,lu):

ave:=subs(x[N]=1,diff(lu,x[N])):

zonot:=[evalf(ave)]:

lu:=lu/x[N]^ave:
lu:=expand(x[N]*diff(x[N]*diff(lu,x[N]),x[N])):
sd:=evalf(sqrt(subs(x[N]=1,lu))):

zonot:=[op(zonot),sd]:

for r1 from 3 to r do
 lu:=expand(x[N]*diff(lu,x[N])):
 zonot:=[op(zonot),evalf(subs(x[N]=1,lu)/sd^r1)]:
od:



lu:=subs({seq(x[i]=1,i=1..N)},gu):
gu:=gu/lu:

for i from 1 to N do
 Ave[i]:=evalf(subs({seq(x[i1]=1,i1=1..N)}
,diff(gu,x[i]))):
od:

gu:=gu/mul(x[i]^Ave[i],i=1..N):

for i from 1 to N do
 SD[i]:=
evalf(sqrt(subs({seq(x[i1]=1,i1=1..N)}, expand(x[i]*diff(x[i]*diff(gu,x[i]),x[i]))))):
od:





menuvalim:=[]:

gu1:=diff(gu,x[N])/SD[N]:

for i from 1 to N-1 do
 kak:=subs({seq(x[i1]=1,i1=1..N)},diff(gu1,x[i]))/SD[i]:
 kak:=evalf(kak):
 menuvalim:=[op(menuvalim),kak]:
od:

zonot,menuvalim:
end:




#StatAnal(L,r):Inputs a list of lists L, all of the same length
#where each entry starts with the name of the individual and
#numerical values of the same attibutes, where the last one is
#the most important.
#Outputs:(1) A list of length nops(L)-1 each of whose entries
#contains [ave,s.d.,alpha_3, ..., alpha_r]
#for the i-th attribute
#(2) A list of lists, let's call it M, such that M[i][j]
#is the correlation between the i-th and j-th attribute
#Try:
#StatAnal(NoHillel(GenData()),10);
StatAnal:=proc(L,r) local gu,x,N,ku,i1,Ave,Mom,lu,Mom1,gu1,SD,i,r1,Cor,j:


if nops({seq(nops(L[i]),i=1..nops(L))})<>1 then
 print(`Bad input`):
 RETURN(FAIL):
fi:

gu:=GF1(L,x):

N:=nops(L[1])-1:

ku:={seq(x[i1]=1,i1=1..N)}:

lu:=subs(ku,gu):
gu:=gu/lu:

for i from 1 to N do
 Ave[i]:=evalf(subs(ku,diff(gu,x[i]))):
od:

gu:=gu/mul(x[i]^Ave[i],i=1..N):

for i from 1 to N do
 SD[i]:=evalf(sqrt(subs(ku,x[i]*diff(x[i]*diff(gu,x[i]),x[i])))):
od:

Mom:=[]:


for i from 1 to N do

 Mom1:=[Ave[i],SD[i]]:
gu1:=gu:
gu1:=x[i]*diff(gu1,x[i]):
gu1:=x[i]*diff(gu1,x[i]):
 for r1 from 3 to r do
 gu1:=x[i]*diff(gu1,x[i]):
 Mom1:=[op(Mom1),evalf(subs(ku,gu1)/SD[i]^r1)]:
 od:
Mom:=[op(Mom),Mom1]:
od:

Mom:

for i from 1 to N do
 for j from 1 to N do
  Cor[i,j]:=evalf(subs(ku,diff(diff(gu,x[i]),x[j]))/SD[i]/SD[j]):
 od:
od:

Mom,[seq([seq(Cor[i,j],j=1..N)],i=1..N)]:


end:




#SipurStatAnalA(L,r,A,Koteret,POPULATION):
#verbose version of StatAnalA(L,r,A)
#where A is a list of attributes (in English)
#and r is a positive integer.
#Try 
#SipurStatAnalA(Etnan(NA),10,
#[Ranking,Tuition,NumberOfStudents,TicketCostForPraying],
#`Statistical Analysis of the Hillel Foundation Ticket Cost`,`US Colleges`);

SipurStatAnalA:=proc(L,r,A,Koteret,POPULATION) local gu,N,r1,i:
gu:=StatAnalA(L,r):
N:=nops(L[1])-1:

print(Koteret):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`We are interested in statistical information about`):
print(A[N]):
print(`in the population`, POPULATION ):
print(`and as it relates to the attributes`):
print(op(1..nops(A)-1,A)):
print(`By collecting data, we found the following data`):
print(`in the format: Name`, op(A)):
for i from 1 to nops(L) do
 print(L[i]):
od:

print(`The average of`, A[N], `over `, POPULATION, `is `):
print(gu[1][1]):
print(`its standard deviation is`):
print(gu[1][2]):

for r1 from 3 to r do
 print(`The `, r1, `-th alpha coefficient is`, gu[1][r1]):
od:


for i from 1 to N-1 do
 print(`The correlation between`, A[N], `and `, A[i] , `is: `):
 print(gu[2][i]):
od:

print(`This ends this fascinating article`):

end:



#SipurStatAnal(L,r,A,Koteret,POPULATION):
#verbose version of StatAnal(L,r,A)
#where A is a list of attributes (in English)
#and r is a positive integer.
#Try 
#SipurStatAnal(NoHillel(GenData()),10,
#[Ranking,Tuition,NumberOfStudents],
#StatisticalAnalysisOfTuitionCosts,UScolleges);

SipurStatAnal:=proc(L,r,A,Koteret,POPULATION) local gu,N,r1,i,j:
gu:=StatAnal(L,r):
N:=nops(L[1])-1:

print(Koteret):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`We are interested in statistical information about`):
print(op(A)):
print(`in the population`, POPULATION ):
print(`and how they relate to each other.`):
print(``):
print(`By collecting data, we found the following data`):
print(`in the format: Name` , op(A)):
for i from 1 to nops(L) do
 print(L[i]):
od:

for i from 1 to N do
print(`The average of`, A[i], `over `, POPULATION, `is `):
print(gu[1][i][1]):
print(`its standard deviation is`):
print(gu[1][i][2]):

for r1 from 3 to r do
 print(`The `, r1, `-th alpha coefficient is`, gu[1][i][r1]):
od:
od:


for i from 1 to N do
  for j from i+1 to N do
    print(`The correlation between`, A[i], `and `, A[j] , `is: `):
     print(gu[2][i][j]):
 od:
od:
print(`This ends this fascinating article`):

end: