Math 640: EXPERIMENTAL MATHEMATICS (Bayesian Probability and Gambling) Spring 2020 (Rutgers University)
http://sites.math.rutgers.edu/~zeilberg/math640_20.html
STARTING March 23, the class in remote via WebEx, per University policy
Last Update: May 20, 2020
 Teacher:
Dr. Doron ZEILBERGER ("Dr. Z")

Classroom:
Allison Road Classroom Building
[Busch Campus], IML Room 119

Time: Mondays and Thursdays , period 3 (12:00noon1:20pm)

[Added March 29, 2020]: Pick a project.

Textbook: "A first course in Bayesian Statistical Methods" by Peter D. Hoff and online resources:

Dr. Zeilberger's Office: HILL 704 [BUT UNTIL FURTHER NOTICE EMPTY]

Dr. Zeilberger's Email: DoronZeil at gmail dot com

Dr. Zeilberger's Office Hours (Spring 2020): Tues. 3:304:30pm; Thurs. 3:304:30pm
DISCONTINUED DUE TO REMOTE LEARNING.
Please use email.
Description
Experimental Mathematics used to be considered an oxymoron, but the future of mathematics is in that direction.
In addition to learning the philosophy and methodology of this budding field, students will become computeralgebra wizards,
and that should be very helpful in whatever mathematical specialty they are doing (or will do) research in.
We will first learn Maple, and how to program with it. This semester we will learn, from an experimental mathematics point of view,
the very important topic of Bayesian Probability, so important in machine learning and data science, as well as
algorithmic gambling.
Optional (but Highly recommended) Getting Started Homework
Teach yourself the basics of Maple by reading Frank Garvan's goldenoldie
part 1,
part 2
Note: a few commands are no longer valid in today's Maple. The
most important one is that " has been replaced by %.
For more details, see this 2002 masterpiece.
How to submit homework

Every class has a homework assignment. Both Mondays and Thursdays' are due
10:00pm Sunday of the following week.
It should be sent to the email address
ShaloshBEkhad at gmail dot com
Subject: HomeWork#X
and then attach a .txt file(s) called
hwXYourFirstNameYourLastName.txt
where X is the number of the assignment
Except for the first assignment, you should have two attached .txt files.
For example, when Yukun Yao submits homework assignments 2 and 3, he should
email ShaloshBEkhad at gmail dot com, with
Subject: HomeWork#2#3
and attach the text files (the source code plus human comments (such lines must start with #)
hw2YukunYao.txt and hw3YukunYao.txt

The very first line of the .txt files should have
#Please do not post homework
or
#OK to post homework
Followed by
#Your Name, Date, Assignment X
Programs done on Thurs., Jan. 23, 2020
C1.txt ,
Contains procedures

MH1(): simulating the MontyHall game without changing

MH2(): simulating the MontyHall game with changing

RS(N): a random sequence of 0 and 1 of length N

FRun(L,k): inputs a list of length L of 0,1, and a pos.integer k, outputs the
location of the end of the first run (if it exists) of k consecutive 0's, or returns FAIL

HotH1(N,k): generates a random sequence of {0,1} of length N, finds
the end location of 0streak of length k, and outputs the value right after it
Homework for Thurs., Jan. 23, 2020, class (due Sunday Jan. 26, 10:00pm)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#1
and an attachment
hw1FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

[People who took this class before, or already know Maple, are exempt from this]
Read and do all the examples, plus make up similar ones,
in the first 30 pages of Frank Garvan's
awesome Maple booklet
(
part 1,
part 2)
Only record a small sample in hw1YourName.txt .

[For everyone].
Read and understand Chapter 1 (the Introduction) of Hoff's book, just to get the flavor of Bayesian thinking.
Do not worry about the details.

[Mandatory for experts, optional for novices]
Simulate the phenomenon of "regression about the mean". Call a sequence of {0,1} of even length N gifted if it has more than N/2+5*sqrt(N) 1's
(i.e. if its sum is more than N/2+5*sqrt(N)). By running RS(N) M times, keep track of the gifted ones, and those that are immediately
followed by another gifted sequence. Convince yourself that most gifted sequences are not followed by gifted ones.
Added Jan. 26, 2020: Read the students' nice solutions
Programs done on Mon., Jan. 27, 2020
C2.txt ,
Contains procedures

Nor1(L): The normalized version of the list of numbers L

LD(L): simulating a loaded die.
Inputs a finite discrete prob. distribution and outputs an integer from 1 to nops(L) (the number of elements in L)
with probability L[i]

SimBay(Pri,Lik,N): simulating Bayes' rule .
Inputs a finite prob. dist., Pri, such that Pr(X=i)=Pri[i] and a list of the same length such that
for some event Pr(SuccessX=i)=Lik[i], finds the APPROXIMATE posterior prob. distribution of 1 through nops(L) if that even happened
by doing it N times.

Bay(Pri,Lik,N): Inputs a finite prob. dist. Pri, such that Pr(X=i)=Pri[i] and a list of the same length such that
for some event Pr(SuccessX=i)=Lik[i], finds the EXACT posterior prob. distribution if that event happened using Bayes' rule

[Added Jan. 28, 2020, thanks to Victoria Chayes]
RSgM(N,L): inputs a positive integer N and a ``loaded die'' L, outputs a random sequence of length N.
Using the Maple Statistics and LinearAlgbera packages.
Homework for Monday, Jan. 27, 2020, class (due Sunday Feb. 2, 10:00pm)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#2
and an attachment
hw2FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

[People who took this class before, or already know Maple, are exempt from this]
Read and do all the examples, plus make up similar ones,
in pages 3060 of Frank Garvan's
awesome Maple booklet
(
part 1,
part 2)
Only record a small sample in hw2YourName.txt .

Read and understand Chapter 2, sections 2.12.4 (pp. 1722) of Hoff's book.

Write a generalization of procedure RS(N) from
C1.txt ,
call it RSg(N,L), that inputs in addition to N, a "loaded die" L.
In particular RSg(N,[1/2,1/2]) should do the same thing as RS(N) (except with 1's and 2's instead of 0's and 1's).

Using the above program RSg(N,L), run the line:
Statistics[Histogram]([seq(add(RSg(3000,L)),i=1..3000)]);
for
L=[1/2,1/2] ( a fair coin), L=[1/3,2/3] ( a loaded coin with Pr(Head)=1/3),
L=[1/6, 1/6, 1/6, 1/6, 1/6, 1/6] ( a fair die),
L=[1/10,1/10,1/10,7/10] ( a loaded tetrahedron)
Look at the resulting shapes. Do they look familiar?
Added Jan. 26, 2020: Victoria Chayes noticed that if you use the simpleminded RSg(N,L) of the previous problem
it would take a very long time. You are welcome to do the above with procedure RSgM(N,L), written by Victoria,
that I added to C2.txt .

[Human exercise] State and prove the fact that if you toss a coin N times with Pr(Heads)=p, the probability that you get exactly k Heads (and hence Nk Tails)
is
N!/(k!(Nk)!)p^{k}(1p)^{Nk}
Added Feb. 2, 2020: Read the students' nice solutions
Programs done on Thurs., Jan. 30, 2020
C3.txt ,
Contains procedures

AgeDS(): The random variable "age" in the people present on Jan. 30, 2020

RSS(N,M): a random probability sample space on {1, ..., N}
with random probabilities (M is a pos. integer)
The output is a list of probabilities such L[i] is the
prob. of the "atomic" event i
M is an arbitrary artificial pos. integer>=1 to generate random prob.

RRV(N,M,K): A random random variable that takes (random) values from K and to K
N and M are as in RSS(N,M)

Av(X):inputs a finite r.v. outputs the expectation (aka average, aka mean)

PGF(X,t): The probability generating function (a finite sum of monomials, or "generalized polynomial"
(t^(3) is OK and t^(1/2) is OK)

Moms(X,k): inputs a finite prob. dist. (with our data structure) and
a pos. integer k, outputs a list of size k whose first entry is
the mean (alias expectation alias average), the second entry is
the variance and the ith entry is the ith moment ABOUT the mean
Homework for Thus., Jan. 30, 2020, class (due Sunday Feb. 2, 10:00pm)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#3
and an attachment
hw3FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

[People who took this class before, or already know Maple, are exempt from this]
Read and do all the examples, plus make up similar ones,
in pages 6090 of Frank Garvan's
awesome Maple booklet
(
part 1,
part 2)
Only record a small sample in hw3YourName.txt .

Read and understand this
this masterpiece

The scaled moments about the mean, of a random variable X is defined by m_r(X)/m_2(X)^(r/2). Use
procedure Moms(X,k) to write a procedure Alpha(X,k) whose first two entries are the same as Moms(X,k)
but the remaining ones are the scaled moments.

We will soon learn that the PGF of the binomial distribution is (p*t+(1p))^{n},
use procedure Moms(X,2) to rigorously find the expectation, and variance,
and the limit of the scaled moments as n goes to infinity of the 3rd through 12th moments.
Added Feb. 2, 2020: Read the students' nice solutions
Programs done on Monday, Feb. 3, 2020
Today was a "traditional class", but Dr. Z. wrote two little procedures
C4.txt ,
Contains procedures

Po(k,a): The Poisson distribution as the limit of B(n,a/n) as n goes to infinity.
k has to be numeric (a positive integer) by a can be symbolic.

cdfBD(n,p,x): The CMF (Comulative Mass function) of the binomial distribution from 0 to trunc(n*p+x*sqrt(n*p*(1p)))
(Recall that the standarddeviation of B(n,p) is sqrt(n*p*(1p))

[Added Feb. 10, 2019, code by Johnny Fonseca, who won $5 dollars (for solving the challenge below]
Po(k,a): The Poisson distribution as the limit of B(n,a/n) as n goes to infinity.
Both k and a are symbolic.
Homework for Monday., Feb. 3, 2020, class (due Sunday Feb. 9, 10:00pm)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#4
and an attachment
hw4FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

[People who took this class before, or already know Maple, are exempt from this]
Read and do all the examples, plus make up similar ones,
in pages 90 to the end of Frank Garvan's
awesome Maple booklet
(
part 1,
part 2)
Only record a small sample in hw4YourName.txt .

The Comulative Mass function of the Binomial Distribution B(n,p) is programmed in
C4.txt ,
as procedure cdfBD(n,p,x): that reads
add(binomial(n,k)*p^k*(1p)^(nk),k=0..trunc(n*p+x*sqrt(n*p*(1p)))):
Investigate how close it is approximated by the CDF of the error function
Int(exp(t^2/2)/sqrt(2*Pi),t=infinity..x)
For p=1/2,x=1, how big does n have to be so the difference of the ratio from 1 is less than 0.001?

The next probability distribution that we will learn about is the Poisson distribution. It
may be defined by procedure Po(k,a) in C4.txt ,
that reads
limit(binomial(n,k)*(a/n)^k*(1a/n)^(nk),n=infinity)
(It is essentially the binomial distribution where the probability of Heads is small, but you toss it many times,
so the expected number of Heads is fixed, a).
This procedure works for symbolic a and numeric k. By running
Po(a,k) for k=1,2,3,4,5,6
Conjecture a pattern for the prob. mass function Pr(X=k) of the Poisson distribution with expectation a.

[Human exercise]: Prove it humanly.

[Optional challenge, I don't know the answer, 5 dollars to the first prover].
Get Maple to prove it by doing Po(a,k) for symbolic a and k.
Added Feb. 19, 2020: Johnny Fonseca won the prize. It is now procedure Po2(a,k) in
C4.txt .
Added Feb. 9, 2020: Read the students' nice solutions
Material covered on Thurs., Feb. 6, 2020
Today we had the honor and pleasure to have a great
guest lecture
by guru Neil Sloane, the founder and president of the
OnLine Encyclopedia of Integer Sequences .
Optional Homework for Thurs., Feb. 6, 2020, class (due indefnitely)

[$50 donation each in honor of the first solvers]
Do assignments HW1 and and HW2 in Dr. Sloane'e lecture
Added Feb. 10, 2020: Look at the students' nice pictures
and partial progress.
Added Feb. 14, 2020: Read Dr. Sloane's message about
Scott Shannon's coloring.
Programs done on Monday, Feb. 10, 2020
C6.txt ,
Contains procedures

Bp(n,p,k): The prob. mass function of the Binomial distribution, tossing a coin n times
Pr(H)=p, k=0..n, Pr(X=k).
Note that n and p are the PARAMETERS and k is the (discrete) `random variable'

HT(n,p,x): Frequentist Hypothesis testing. If the owner of the coin claims that Pr(H) is at least
p what is the cutoff of number of Heads below which you can be confident with level of confidence
x that he is a lier.

dbeta(p,a,b): The Beta distribution with parameters a and b.
Note that a and b are the PARAMETERS and p is the (continuous) `random variable'
Homework for Monday., Feb. 10, 2020, class (due Sunday Feb. 16, 10:00pm)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#6
and an attachment
hw6FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

Write a procedure:
CheckHT(n,p,x,N), that first calls HT(n,p,x) to find out what is the "cutoff" below which we have to reject the "null hypothesis" that
the coin has Pr(Heads)=p with (onesided) confidence level x, and then runs the line
with(Statistics): convert(Sample(RandomVariable(BernoulliDistribution(p)),n), + );
N times (to find out the number of Heads) in n random coin tosses with Pr(Heads)=p), and computes the ratio of those for
which you have to reject the Null Hypothesis. Find out how close there are. Experiment with p=1/6,1/3,1/2,2/3 and n=600,
and x=0.95, x=0.99, and N=1000:

Consult the
Maple Statistics help page
to find out how to access the Binomial distributions.
Run it many times and examine the histograms and compare it to the theoretical histogram according to binomial(n,k)*p^k*(1p)^k

Read section 3.1 of Hoff's book

Write a procedure
BCI(f,t,x)
(Bayesian confidence interval)
that inputs a continuous pdf f , in the variable t defined on [0,1] (e.g. the Beta distributions)
and a positive integer x (in application close to 1) and outputs and outputs a symmetric interval,
containing 1/2, of the form [1/2z,1/2+z], that contains a portion x of the distribution, i.e. such that
int(f,t=1/2z..1/2+z)=x
What are
BCI(dbeta(a,a,t),t,x) for a=1,2,3,4, x=0.95, 0.99, 0.999)

[Challenge, I don't know the answer, maybe it does not exist]
Find a formula in terms of a and x for BCI(dbeta(a,a,t),t,x)

[Challenge] Write a procedure
HPD(f,t,x)
that inputs f, t, and x as above but outputs a HPD interval (see Def. 7, bottom of p. 42).
Added Feb. 16, 2020: Read the students' nice solutions
Programs done on Thurs., Feb. 13, 2020
C7.txt ,
Contains procedures

PostD(Pri,t,a,b,Lik,k,Data1): inputs a CONT. pdf Pri in the variable t,
for the Prior distribution, whose domain is t=a to t=b, a
DISCRETE r.v. Lik in terms of k, and a LIST of the data, Data1, uses
the Bayesian methodology to output the PDF of the POSTERIOR distribution
For example,
PostD(t^(31)*(1t)^(31),t,0,1,t^k*(1t)^(1k),k,[0,1,0,1,0]);

dgamma(t,a,b): The Gamma distribution
Khizar Anjum has prepared , after this class, this
useful Maple code .
Special Homework for Thurs., Feb. 13, 2020, class (due Friday Feb. 14, 10:00pm)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#VD
and attachments
hwVDFirstNameLastName.txt
and
hwVDFirstNameLastName.jpg
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

[A chocolate box to the best entry voted by students]
(i)Write a procedure
IsGood(c)
that inputs a complex number and decides whether iterating the
Mandelbrot map z>z^2+c, stating at z=0 is bounded or not.
(ii) By taking a fine enough grid in the complex plane, 3 ≤ x, y ≤ 3, with resolution 1/M write a procedure
GoodPoints(M)
that outputs the good points
(iii) By using plot(SetOfPoints, style=point) draw (an approximation to) the famous Mandelbrot set
(iv): It looks a little like a sideway Valentine. Rotate to make it look like a Valentine.
(v) Make it into a .jpg file and send me the attachment AND the code.
Added Feb. 16, 2020: Read the students' nice pictures and Maple code.
Robert DoughertyBliss won the chocolate box for this amazing Zalentine.
Regular Homework for Thursday, Feb. 13, 2020, class (due Sunday Feb. 16, 10:00pm)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#7
and an attachment
hw7FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

Read Chapter 4 of Hoff's book

Implement the fourth line from the bottom on p. 54 of Hoff's book, by writing a procedure
that approximates integrals with respect to a known pdf for which Maple knows to sample,
like I tried to do in class.

Explain what I did wrong.
Added Feb. 13, 2020, 9:30pm: I figured it out, I think. To get an approximation to the integral of
Int(exp(x^2/2)/sqrt(2*Pi)*f(x),x=infinity..infinity)
You simply do, if L=Sample(RandomVariable(Normal(0,1)),N);
add(f(L[i]),i=1..N)/N:
Added Feb. 16, 2020: Read the students' nice solutions
Stuff done on Monday, Feb. 17, 2020
Today started with voting on the Valentine day competition,
won by Robert DoughertyBliss .
Then we explored some of the builtin random variables in the Maple Statistics package, and compared, exact integration,
numerical integration (using evalf(Int)) and approximating the integrals by Monte Carlo.
Homework for Monday, Feb. 17, 2020, class (due Sunday Feb. 23, 10:00pm)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#8
and an attachment
hw8FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

Write a procedure
MyMonte(f,t,a,b,N)
that inputs a function f of t (given as an expression in t), that has support in the interval a ≤ t ≤ b (i.e. is 0 otherwise),
and that is always positive in a < t < b, and whose integral from a to b is 1 (in other words, it is a pdf),
and outputs N random samples from the random variable whose pdf is f. Use the method at the bottom of p. 385
Lambert's great book.
If the CDF is not closedform, or its inverse function of the CDF is not closed form it should return FAIL.

Use the above procedure, MyMonte(f,t,a,b,N), to write a procedure
MonteInt(f,g,t,a,b,N)
that uses MyMonte(f,t,a,b,N) with N big to approximate the integral
Int(g(t)*f(t),t=a..b) . Compare the output to the exact value (when possible).
What are
(i) MonteInt(2*t,sqrt(1+t),t,0,1,N)
(ii) MonteInt(2*t,sin(t),t,0,1,N)
(iii) MonteInt(5*t^4,cos(t^3+1),t,0,1,N)

Use Maple's Statistics's commands RandomVariable(Normal(a,b)) and Sample, to write a procedure
F(a1,b1,a2,b2,f,g,x,y,N)
That inputs real numbers a1,b1,a2,b2, and an expression f in x,and an expression g in x and y,
and a large positive integer N, and outputs an approximation to the
double integral
1/(2*Pi)*Int(g*exp((xa1)^2/(2*b1)*exp((ya2)^2/(2*b2),R)
where R is the region
{infinity < x,y < infinity ; f(y) < x} .
For random values of the parameters (including f and g) compare the output to the exact result.

Write a procedure
AvMax(n,N)
that uses N random samples of NormalDistribution(0,1) (the standard normal distribution), and outputs an approximate for
the expectation of the r.v.
max(X1,X2, ..., Xn)
check that the values of AvMax(n,10^7) for n=2,3,4,5 are close to those listed in Marcus Michelin's
nice note (last page).
Try to compute approximations to max(X1, .., Xn) for n from 6 to 10 to three decimal digits, by making N large enough.

[Human/machine Challenge, $50 donation to the OEIS in honor of the first discoverer] Find the exact value, in terms of Pi, e, and any of the standard mathematics function for the case n=6.

[Human/machine Challenge, $150 donation to the OEIS in honor of the first discoverer] Find an efficient algorithm to find the exact value for
max(X1, ..., Xn) for any n.
Added Feb. 23, 2020: Read the students' nice solutions
Programs done on Thurs., Feb. 20, 2020
C9.txt ,
Contains procedures

RRV(K1,K2a,K2b,N): inputs pos. integer K1, and integers K2a, K2b, and a large
pos. integers N, and outputs a random discrete random variable of size N with the format
[[value of X, prob. of ], ...]

EAv(X,f,x): inputs a discrete r.v. X, and an expression f in terms of the variable x.
Outputs the EXACT value of E[f(X)].

DAv(X,f,x,M): An estimate for EAv(X,f,x) (see above) by
Democratic sampling (not "importance" sampling) using M random
places (regardless of their "impotance")

NickAv(X,f,x,M1,M2): An estimate for EAv(X,f,x) (see above) by
Importance sampling, using the famous Metropolis algorithm using M1 "warm up"
rounds (w/o taking records) followed by M2 real rounds
places (paying full attention to of their importance).
Corrected after class, thanks to Jingyang Deng

NickAvWrong(X,f,x,M1,M2): Original WRONG wrong version of NickAv(X,f,x,M1,M2) with the inequality r < rat reversed
[Only to be used in homework below]
Homework for Thurs. Feb. 20, 2020, class (due Sunday Feb. 23, 10:00pm)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#9
and an attachment
hw9FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

[Challenge, I don't know the answer]: In class, I entered the inequality deciding whether to reject or accept the "next guy".
in the wrong direction. It is preserved in procedure
NickAvWrong(X,f,x,M1,M2):
Surprisingly, we still got pretty good approximations with the wrong version. Run it many times, for many random choices of
RRV(K1,K2am,K2b)
and compare the exact value of the expecation, the "democratic" version (DAv) , the correct Metropolis, NickAv (correct by Jingyang Deng),
and the wrong version NickAvWrong. Is NickAv significantly better than NickAvWrong?

Wtite a continuous analog of EAv(X,f,x), call it CAv(P,t,f,a,b)
that inputs a NONnormalized probability density function P given as an expression in t,
another expression in t, called f, and real numbers a and b (possibly infinite)
such that P and f are supported in the interval a < t < and that computes (exactly, either with int, or using numerical integration with Int and evalf)
the integral of P*f from t=a to t=b divided by the integral of P from t=a to t=b

[A bit of a challenge] Write a continuous analog of NickAv, call it
CNickAv(P,t,f,a,b,M1,M2)
that uses the Metropolis algorithm with a warm up period of M1 rounds, and a "datacollecting" period of length M2, to
esimate the exact value of CAv(P,t,f,a,b).

What do you get for example
CNickAv(1/(t^3+1),t,t,0,10,1000,10000) vs CAv(1/(t^3+1),t,t,0,10)
CNickAv(exp(t^3),t,t^2,0,10,1000,10000) vs CAv(exp(t^3),t,t^2,0,10)
Experiment with other random choices.
Added Feb. 23, 2020: Read the students' nice solutions
Programs done on Monday, Feb. 24, 2020
C10.txt
Contains procedures:

IterS(f,x,x0,eps,N): inputs an expression f in x, an initial value x0, a small positive number, eps, and a large positive integer N,
outputs an approximation (within eps) to the equation
f(x)=x
By starting at x=x0 and iterating N times, or returns FAIL.
For example, to approximate sqrt(2)1 try:
IterS((1x^2)/2,x,0.4,1000)

MIterS:(F,x,x0,eps,N): Solving a SYSTEM var=F(var) where F is a list of functions
and x is a list of variables,e.g. x=[x1,x2]. Try:
MIterS([(1+x+y)/(2+x+y), (5+x+y)/(2+3*x+y)],[x,y],[1,1],10^(7),100);

GenBiPoisson(L1,L2,n,N): inputs positive integer N, another pos. integer
n, such that 1 ≤ n ≤ N, and two positive real numbers L1, L2, outputs a list
of length N whose first n entries are sampled from Poisson(L1), and
the remaining ones from Poisson(L2)
Reading homework for Monday Feb. 24, 2020, class (due Thurs. Feb. 27, 2020)

Read and understand
Bayesian Inference: Gibbs Sampling
by Ilker Yildirim
Homework for Monday Feb. 24, 2020, class (due Sunday March 1, 10:00pm)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#10
and an attachment
hw10FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

[A human challenge] : Procedure IterS(f,x,x0,eps,N), sometimes converges to a root of the equation f(x)=x if x0 is close enough to it
and sometimes it does not. Can you predict, theoretically, if the root is r, whether it does or does not?

[A bigger human challenge] Procedure MIterS(F,x,x0,eps,N), sometimes converges to a root of the system F(x)=x, let's call it r (r is a list of numbers)
if x0 is close enough to x0. Can you predict, a priori, when it does?
[Hint: Jacobian, eigenvalues]

After class I noticed that the code of
IterS(f,x,x0,eps,N)
does not reflect Gibbs sampling accurately, as described in Professor Yildirim's
Bayesian Inference: Gibbs Sampling
(and also in Hoff's book). There you don't update the whole list at once, but bulidup, using a combination of
the old values and the previous new values.
Write a procedure
BetterIterS(f,x,x0,eps,N)
that does it that way.
Added March 2, 2020: Read the students' nice solutions
Programs done on Thurs., Feb. 27, 2020
C11.txt
Contains procedures:

GD1(p1,p2,n,N): "Synthetic data" (for testing) Bernoulli(p1) for the first N tosses
followed by Bernoulli(p2) for the last Nn tosses, with ONE game

GD(p1,p2,n,N,M): generates M samples from GD1(p1,p2,n,N),

UnNorPost(DA,p1,p2,n): The unormalized posterior trivariate joint prob. distribution (p1,p2,n) deduced from the data DA
assuming uniform priors for both p1 and p2, and uniform prior for n. The output is symbolic in p1,p2, but
the dependence on the discrete (finite) variable n is given as a list of length N (N=nops(DA[1]))

Post(DA,p1,p2): The posterior distribution from the data DA, given
as list of length N, whose nth item is the pdf for p1,p2 with crossover n
asuming uniform priors for both p1 and p2, as well as n.

MARGp1(DA,p1): Inputs a data set DA, and a symbol p1, outputs the EXACT (symbolic)
expression in p1 for the MARGINAL posterior of p1 (where p2 is integrated over, and n is summed over)
followed by the Expected (predicted) value for p1

MARGp2(DA,p2): Inputs a data set DA, and a symbol p12, outputs the EXACT (symbolic)
expression in p2 for the MARGINAL posterior of p2 (where p1 is integrated over, and n is summed over)
followed by the Expected (predicted) value for p2

MARGn(DA): Inputs a data set DA,, outputs the probability mass function supported between 1 and N ( that nops(DA[1])),
given as a list of length N, of the Marginal n
(where p1 and p2 are "averaged out", i,e. intergated over), followed by the expected (i.e. predicted) value of the crossover point,
judging from the data
Homework for Thurs. Feb. 27, 2020, class (due Sunday March 1, 10:00pm)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#11
and an attachment
hw11FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

Implement the NONSAMPLING part of
Bayesian Inference: Gibbs Sampling
by Ilker Yildirim,
by writing "Poisson Analogs" for all the procedures in
C11.txt
Now the parameters are Lambda1, Lambda2 (call them L1,L2), rather than p1 and p2, as well as n.

Extend the programs in
C11.txt
to the situation where there are two crossovers, n1,n2, where the bias of the coin from 1 to n1 is p1, from n1+1 to n2 is p2,
and from n2+1 to N is p3. You need MARGp1, MARGp2, MARGp3, MARGn1, MARGn2, and all the sipporting procedures.
For some random data sets, plot the marginal distributions.

[Optional Challenge, a prize for the best entry]
Implement Gibbs Samplig for the Bernoulli case we did in class. The sampling distribution for p1 and p2 are
(conditioned on n) the appropriate Beta Distribution, but the sampling of n (conditioned in p1 and p2)
is a PROBABILITY table, that may cause trouble the way Maple handles floating points (and I was unable
to do it for the Poisson case). Find approximations to the expected posterior p1, p2, and n and compare them
to the "exact" values given by Bayes.
[Added Feb. 29, 2020:
Note that the output from
Sample(RandomVariable(ProbabilityTable(L)),1)[1] is "floatingpoint", so you need to convert it into
an integer by doing trunc().]

[Bigger optional Challenge, a prize for the best entry]: Do Gibbs sampling for the Poisson case.
Here the probability table is inherently floatingpoint, so you need to convert it into rational
numbers, and then normalize so that you the probabilities addup to 1.

[General, VERY OPTIONAL, BIG challenge, unrelated to the class, but looks interesting. No deadline] Study this
problem,
by Rutgers PhD alumnus Hanlong Fang.
Added March 2, 2020: Read the students' nice solutions
Programs done on Monday, March 2, 2020
C12.txt
Contains , in addition to all the procedures from C11.txt, (use Help11() so see them), the following procedures:

ParamBetaP1(DA,n,alpha,beta): The paramaters [alpha',beta'] of the Beta distribution for
p1 if the the prior is Beta(alpha,beta), the data set is DA and the crossover is n

ParamBetaP2(DA,n,alpha,beta): The paramaters [alpha',beta'] of the Beta distribution for
p2 if the the prior is Beta(alpha,beta), the data set is DA and the crossover is n

PTn(DA,p1,p2): The probability table of length N that gives you the prob. dist.
of the crossover point being i, assuming that you know p1, and p2

GibbsSamp(DA,M1,M2,alpha,beta): Gibbs Sampling for the crossover problem with
Bernoulli trials, with burnout period M1, followed by M2 "for real" samples

AvBay(DA): The predicted average of n, p1, and p2, from the data set DA with alpha=beta=1
ACCORDING TO BAYES ("exactly")

AvGibbs(DA,M1,M2): The predicted average of n, p1, and p2, from the data set DA with alpha=beta=1
using Gibbs Sampling with burnout M1 and M2 samples
Homework for Monday, March 2, 2020, class (due Sunday March 8, 10:00pm)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#12
and an attachment
hw12FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

Modify procedure AvBay(DA) to write a procedure
MomBay(DA,k)
that inputs a data set DA, as before, and a positive integer k, and outputs a list of length 3 whose
first entry is the list of length k whose first entry is the expected n (as before), the second entry is its variance, the the jth
entry is the jth entry about the mean of n
Second entry is the list of length k whose first entry is the expected p1 (as before), the second entry is its variance, the the jth
entry is the jth entry about the mean of p1
Third entry is the list of length k whose first entry is the expected p2 (as before), the second entry is its variance, the jth
entry is the jth entry about the mean of p2
[Hint: For a continuous random variable with pdf f(x), the average is mu=Int(x*f(x),x), the variance is Int((xmu)^2*f(x),x), and more
generally, the kth moment about the mean is Int((xmu)^k*f(x),x). Analogously for the discrete case]

Modify procedure AvGibbs(DA,M1,M2) to write a procedure
MomGibbs(DA,M1,M2,k)
that estimates the output of MomBay(DA,k). For random choices of DA (using
GD(p1,p2,n,N,M)), compare the outcomes of MomBay(DA,k) and MomGibbs(DA,M1,M2,k)
(with different M1's and M2's)

Do the Poisson analog of what was done in C12.txt ,
i.e. Gibbs Sampling for the crossover point, where instead of Benoulli you have Poisson, and instead of
of the Beta distribution you have the Gamma distribution.
Added March 9, 2020: Read the students' nice solutions
Programs done on Thurs., March 5, 2020
C13.txt
Contains ,

RM(M,N): a random M by N matrix whose entries are taken from {1,1} given as a list of lists

PlotMat(A): plots the matrix A where 1 is black and 1 is red

Ene(A): The energy of the matrix A of {1,1}

Mag(A): the magnetization of the spinconfig. A

Vecs(N): the set of all 2^N {1,1} lists of length N

Mats(M,N): the set of all M by N {1,1} matrices

ParF(M,N,x,y): the Gibbs partition function of the M by N Ising model, where x and y are physical variables related to the
absolute temperature and EXTERNAL magenetic field
Homework for Thurs., March 5, 2020, class (due Sunday March 8, 10:00pm)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#13
and an attachment
hw13FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

Write a procedure
RMg(M,N,S)
that inputs pos. integers M and N and a finite set S, and outputs a random M by N matrix whose entries are from S

Write a procedure
PlotMatG(A,S): that inputs a matrix A (given as a list of lists) and a set S of pairs {[s,color]}, where the matrix A
has entries from the set of first compoents of S and s is replaced by its color.

[A bit of challenge]
Write a procedure
ArraysS(S,L)
that inputs a set S and a list L of pos. integers and outputs the set of
L[1] x L[2] x ... x L[nops(L)]
The elements of the output should be listsoflistsoflistsoflists, ... nested k times. for example one member
of ArrayS({0,1,2},[2,2,2]) is [ [[0,2],[1,2] ], [ [1,1],[2,2] ] ].

Read again
Chapter about the Metropolis Algorithm
in Julia Yeomans great book "Statistical Mechanics of Phase Transitions.
and, if possible,
Kenneth G. Wilson's seminal article
(Sci. Amer. 1979)
Added March 9, 2020: Read the students' nice solutions
Programs done on Monday., March 9, 2020
C14.txt
Contains ,

Neis(c,M,N): Inputs a location c=[i,j] and pos. integers M and N outputs the SET
of four neighbors of the cell [i,j], {[i \pm 1,j ], [i, j \pm 1] but with toral convention
[CODE by Alison Bu and Yukun Yao, with improvements by Charles Kenney]

NewA(A,c): The result of flipping the spin at location c in the matrix A
followed by the change in Energy, and change in Magnetization

OSM(A,M,N,x,y): one step in the Metropolis algorithm inputs an M by N matrix of {1,1}
and pos. numbers x,y outputs one step in the Metropolis algorithm according to
the importance x^Ene(A)*y^Mag(A)

Nick(M,N,x,y,K1): Starts with a random M by N {1,1} matrix with physical parameters x and y
and a positive integer K1 runs Metropolis K1 times
Homework for Monday, March 9, 2020, class (due Sunday March 22, 10:00pm, but it is recommended to do it by March 15)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#14
and an attachment
hw14FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

Write a procedure
Neis3D(c,M,N,K) that inputs a location in the form [i,j,k] and pos. integers M,N,K, and outputs the set
of its eight neighbors in a "wrappedaround" M by N by K box.

Wtite a procedure
NewA3d(A,c)
that inputs a threedimensional array whose entries are {1,1} (i.e., a 3D configuration of spins) and a location c in
the form [i,j,k], and outputs the matrix obtained by filpping the spin in that location. It should also return the
change in energy, and the change in magnetization

Write procedures OSM3d(A,M,N,K,x,y) and Nick3D(M,N,K,x,y,K1) that are 3D analogs of procedures OSM(A,M,N,x,y,K1) and Nick3D(M,N,x,y,K1)
Added March 22, 2020: Read the students' nice solutions
Programs done on the reomote class of Thurs., March 12, 2020
C15.txt
[Corrected March 15, 2020]
Contains ,

ExactAv(M,N,x,y,G,En,Ma): Inputs positive integers M and N and real numbers x and y, positive integers K1 and K2
and an expression G that is a MONOMIAL in the symbols E and M (denoting energy and magnetization, respectively), outputs the
the EXACT average of G(En,Ma) over all 2^(M*N) possible M by N {1,1} matrices with the Ising prob. distribution with
temperature parameter x and external field parameter y. It uses "generatingfunctionology". If the
Prob. gen. fun of a r.v. is f(z) then its expectation is (z*d/dz f(z))z=1 and its rth moment is
(zd/z)^r) f(z) z=1 , similarly for more than variable.
It uses generating functions.
For example, to the exact average of the energy with x=1.1 and y=1 on the 4 by 4 rectangular torus do:
ExactAv(4,4, 1.1, 1.0, En,En,Ma).
To get the exact average of En^2*Ma^2 do
ExactAv(4,4, 1.1, 1.0, En^2*Ma^2,En,Ma).

ExactAvBrute(M,N,x,y,G,En,Ma): Same as ExactAv(M,N,x,y,G,En,Ma) but by brute force.

NickAv(M,N,x,y,K1,K2,G,En,Ma):
Inputs positive integers M and N and real numbers x and y, positive integers K1 and K2
and an expression G in the symbols E and M (denoting energy and magnetization, respectively), outputs the
estimate produced by running the Metropolis algorithm K1 in the "warmup" stage (w/o taking records)
followed by K2 iterations where a record of G(E,M) is taken, and added to the average.
For example, to get an estimate for the average of the energy, (where the numerical value if the physical parameters x and y are 2.0 and 1.1) try:
NickAv(4,4, 2.0 ,1.1, 100,1000, E, E,M);
and to get an estimate for the average of the energysquared times the magnetization, type
NickAv(4,4,2.0 ,1.1, 100,1000, E^2*M, E,M);
Homework for Thurs. March 12, 2020, class (due Sunday March 22, 10:00pm, but it is recommended to do it by March 15)
Version of March 13, 2020: (The previous code had an error, and one of the previous questions was deleted)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#15
and an attachment
hw15FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

Read and understand the three new procedures in
C15.txt , and
experiment with them. Try to do:
NickAv(4,4, 1.3 , 1.0, 1000,10000, E, E,M);
ten times, and see whether they are close (when I tried it, it was a little disappointing)

Now try:
NickAv(40,40, 1.3 , 1.0, 1000,10000, E, E,M);
ten times. Are the outputs closer to each other?

According to Julia Yeomans's book starting with either the all 1matrix or the alternationg 1,1, rather than
a random matrix gives better results. Write procedures
NickAvONE(M,N,x,y,K1,K2,G,En,Ma):
NickAvALT(M,N,x,y,K1,K2,G,En,Ma):
that instead of starting with a random matrix, starts with the all 1s and the 1,1,... matrices respectively,
and see if you get better agreements between the various runs.

Write a procdeure
Shrink(A)
that starts with a 3^{k} by 3^{k}
matrix of {1,1} and outputs a 3^{k1} by 3^{k1} matrix, that returns the
"majority rule" in each 3 by 3 consecutive square.
Run
Nick(81,,81,x,y,10000):
for various numerical values of x (not too far from 1), then shrink them successively, and use
PlotMat(A) to plot them. Try to reproduce the kind of pictures in
Kenneth G. Wilson's seminal article
Added March 22, 2020: Read the students' nice solutions.
See also the beautiful
pictures.
Programs discussed during the reomote class of Monday, March 23, 2020
C16.txt
Contains

StatAnal(P,t,K): Inputs a polynomial (or even a function whose Taylor expansion has positive coefficients)
describing the (notnecessarily normalized) probability generating function of some discrete random variable
outputs a list of length K whose
(i) First entry is the expectation (alias mean, alias average)
(ii) Second entry is the variance (aka as second moment about the mean)
(iii) thirdthrough K entries are the scaled moments about the mean
in particular the third entry is the skewness and the fourth moment is the kurtosis. For example, try:
StatAnal((1+t)^1000,t,6);

GRsim1(N,p,n): Simulating ONE Gambler's ruin game with initial capital of n dollars,
At each turn the gambler's win a dollar with probability p and loses a dollar with probability 1p
the game ends as soon as it reaches N dollars or 0 dollars (when the gambler's is ruined).
The output is a pair [0,m] or [1,m], according to whether the gambler's left the casino broke or
rich, and m is the number of rounds. Try:
GRsim1(10,1/2,5);
GRsim1(10,10/21,5);

GRsim(N,p,n,M,s,t): Inputs a positive integer N, a rational number p between 0 and 1, and a pos. integer n<=N
as well as symbols (variables) s and t, outputs
The Empirical probability generating function of running GR1(N,p,N) M times
it returns a polynomial of degree 1 in s, i.e. of the form F(t,s)=P0(t)+s*P1(t) where
PO(t) records the losing cases and P1(t) the winning cases according to duration.
In particular the coefficient of s in F(1,s) is the prob. of winning
the coefficient of s^0 (the constant term) in F(1,s) is the prob. of losing
and F(1,t) is the prob. generating function for the duration of the game (regardless whether you won or lost)
Try:
GRsim(40,1/2,20,100,s,t);

[Added after class, code by Victoria Chayes]
PH(N,p,n): Drawing ONE Gambler's ruin game with initial capital of n dollars. Try:
PH(10,1/2,5);
PH(10,10/21,5);

[Added after class]:
Roulette(L): Inputs a list of pairs [Gain,ProbOfGain], returns Gain with probability ProbOfGain.
For example, for the classical Gambler's ruin problem with probability of winning a dollar being p
and losing a dollar being 1p, you do Roulette([[p,1], [1p,1]]);
If you want to win 6 dollars with prob. 1/2, win 3 dollars with prob. 1/3 and lose 20 dollars with prob. 1/6.
you do
Roulette([[1/2,6],[1/3,3],[1/6,20]]);
Homework for Monday, March 23, 2020, class (due Sunday March 29, 10:00pm)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#16
and an attachment
hw16FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

Read, understand, and experiment with, the two procedures (the first one by Victoria Chayes), PH(N,p,n)
and Roulette(L).

To get statistical analysis about the duration of Gambler's ruin game with prob. of winning a dollar
being p, and with initial capital n, running it M times, regardless of whether you won or lost,
you do
StatAnal(GRsim(N,p,n,M,1,t),t,K)
where K is the number of moments you are interested in.
Run
StatAnal(GRsim(30,p,15,5000,1,t),t,8);
for p=1/10,2/10,3/10,4/10,5/10, and observe the values, especially of the expectation.

Generalize GRsim1(N,p,n) and GRsim(N,p,n,M,s,t) to write, procedures
GenGRsim1(N,L,n) and GenGRsim(N,L,n,M,s,t)
where L is the roulette.
[Added March 26: The convention is that the first time you have a nonpositive amount or
an amount >=N then the game is over
[Note that GenGRsim1(N,[[p,1],[1p,1]],n) is equivalent to GR1sim(N,p,n)]

Run
StatAnal(GenGRsim(N,L,n,10000,1,t),t,6);
For

L=[[1/3,3],[1/6,6],[1/2,4]]

L=[[1/3,3],[1/6,6],[1/2,3]]

L=[[1/3,3],[1/6,6],[1/2,5]]
Do you notice considerable difference between the first and the second outcome?
Between the second and the third outcome? Why?

Convince yourself, that to get the approximate probability of winning the (generalized) game with a roulette, L,
by simulating it M times
you do
coeff(GenGRsim(N,L,n,M,s,1),s,1):
Find these (approximate) probabilities for the above three roulettes.
Do you notice analogous differences?
Added March 29, 2020: Read the students' nice solutions.
Programs discussed during the reomote class of Thurs., March 26, 2020
C17.txt
Contains

ProbWinning(N,p): inputs a positive integer N and a prob. p (either numeric or symbolic),
outputs a list of length N1 whose ith entry is the probability of exiting a winner in
Gambler's ruin game where you start with a capital of i dollars, and then at each round
you win a dollar with prob. p and lose a dollar with prob. 1p. Try:
ProbWinning(10,1/2);
ProbWinning(10,1/3);
ProbWinning(10,p);

ExpDuration(N,p): inputs a positive integer N and a prob. p (either numeric or symbolic),
outputs a list of length N1 whose ith entry is the expected duration of the a
Gambler's ruin game where you start with a capital of i dollars, and then at each round
you win a dollar with prob. p and lose a dollar with prob. 1p. Try:
ExpDuration(10,1/2);
ExpDuration(10,2/3);
ExpDuration(10,p);

PGF(N,p,t): inputs a positive integer N and a prob. p (either numeric or symbolic),
and a variable name t,
outputs a list of length N1 whose ith entry is the PROBABILITY GENERATING FUNCTION FOR DURATION
(the rational function in t whose coefficient of t^i is the probability that the game will end exactly after i round)
In a Gambler's ruin game where you start with a capital of i dollars, and then at each round
you win a dollar with prob. p and lose a dollar with prob. 1p. Try:
PGF(10,1/2,t);
PGF(10,1/3,t);
PGF(10,p,t);

ProbWinningCF(N,n,p): inputs SYMBOLS n, N and a prob. p (either numeric or symbolic),
outputs the CLOSEDFORM EXPRESSION for the prob. of winning in
Gambler's ruin game where you start with a capital of n dollars, and then at each round
you win a dollar with prob. p and lose a dollar with prob. 1p. Try:
ProbWinningCF(N,n,p);

ExpDurationCF(N,n,p): inputs SYMBOLS n, N and a prob. p (either numeric or symbolic),
outputs the CLOSEDFORM EXPRESSION for the Expected duration
Gambler's ruin game where you start with a capital of n dollars, and then at each round
you win a dollar with prob. p and lose a dollar with prob. 1p. Try:
ExpDurationCF(N,n,p);

PGCcf(N,n,p,t): inputs SYMBOLS n, N and a prob. p (either numeric or symbolic), and a variable name (symbol) t;
outputs the CLOSEDFORM EXPRESSION for the Prob. Generating function in t, of the duration
Gambler's ruin game where you start with a capital of n dollars, and then at each round
you win a dollar with prob. p and lose a dollar with prob. 1p. Try:
PGFcf(N,n,p,t);
Homework for Thurs., March 26, 2020, class (due Sunday March 29, 10:00pm)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#17
and an attachment
hw17FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

Using procedure StatAnal(P,t,K) of
C17.txt
(not C16.txt), because I modified it replacing limit(..., t=1) instead of subs(t=1,...) that had issues of dividing by 0,
find closed form expressions in terms of the symbols k and N for the expectation, variance, skewness and
Kurosis for the duration of a Gambler's ruin problem with a fair coin starting with k*N dollars and
exiting when you either have 0 dollars or N dollars. What is the limit of the variance dividied by N^2 as N goes to infinity?
What is the limit (if it exists) in terms of k of the skewness (scaled thirdmoment about the mean),
the kurtosis (scaled fourthmoment about the mean)?

[A bit of a challenge]
Generalize procedure PGCcf(N,n,p,t), call it
PGCcfGen(N,n,p,t,K), that also inputs a positive integer K and finds
the prob. gen. function for the duration of a generalized Gambler's ruin problem
where
You win i dollars with probability p[i], for K<=i<=K, where (it is assumed that they addup to 1, but Maple does
not have to be notified). The game ends as soon as the capical is nonpositive or strictly larger than N1.
[Of course "winning" a negative amount is losing].
Check that your output for K=1, and p[1]=1p, p[0]=0, p[1]=p agrees with the output of GCFcf(N,n,p,t)
[Corrected March 30, 2020 (Thanks to Johnny Fonseca]
Hint: You still use rsolve, but now you need more initial conditions, that later on you have to solve for using
the boundary conditions at the other end. The order of the recurrence is 2*K, the initial conditions
are x[0]=1, x[1]=1, ... x[(K1)]=1 and the boundary conditions are
x[N]=1, x[N+1]=1, ..., x[N+K1]=1. (Since t^{0}=1)
Added March 29, 2020: Read the students' nice solutions.
Programs discussed during the reomote class of Mon., March 30, 2020
C18.txt
Contains

StatAnalPureL(P,t,K): The expecation, variance, and the UNSCALDED moments about the mean
the first entry is the expectation, the second the variance, and then the third, ..., up to the
Kth moment about the mean.
It uses limit(...,t=1) instead of subs(t=1, ...) MUCH SLOWER
Try:
StatAnalPureL((1+t)^n,t,10);

StatAnalPure(P,t,K): The expecation, variance, and the UNSCALDED moments about the mean
the first entry is the expectation, the second the variance, and then the third, ..., up to the
Kth moment about the mean. Try:
StatAnalPure((1+t)^n,t,10);

GP1(L,n,d): inputs a list L , a variable name n, and a positive integer d, guesses a polynomial of degree ≤ d in n,
let's call it P(n) such that P(i)=L[i]. The length of L must be at least d+3. If it fails it returns FAIL. Try:
GP1([seq(i^3,i=1..10)],n,3);

GP(L,n): inputs a list L , a variable name n, guesses a polynomial of degree <=nops(L)3 in n,
let's call it P(n) such that P(i)=L[i] for all i. If it fails it returns FAIL. Try:
GP([seq(i^3,i=1..10)],n);

MomsGR1(n,N,K): inputs a variable name n, and a positive integer N (sufficiently large) and a positive integer K
tries to find polynomial expression in n for the first K moments of the r.v. duration of a Gambler's
Ruin game with a fair coin with exit capital N and starting capital n. If it can only partially do it,
it returns what it can. Otherwise it returns all the full list of K moments. Try:
MomsGR1(n,20,6);

MomsGR(n,N,K,N1,N2): tries to conjecture polynomial expressions in n and N for the first K moments,
by "colleting data from N=N1 to N=N2, where N1 and N2 are positive integers. Try:
MomsGR(n,N,4,12,24);
Homework for Monday, March 30, 2020, class (due Sunday April 5, 10:00pm)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#18
and an attachment
hw18FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

[Recycled from the last homework. Few people got it right, mainly because it was badly stated and
the hint had an error]
Generalize procedure PGCcf(N,n,p,t), call it
PGCcfGen(N,n,p,t,K), that also inputs a positive integer K and finds
the prob. gen. function for the duration of a generalized Gambler's ruin problem
where
You win i dollars with probability p[i], for K<=i<=K, where (it is assumed that they addup to 1, but Maple does
not have to be notified). The game ends as soon as the capical is nonpositive or strictly larger than N1.
[Of course "winning" a negative amount is losing].
Check that your output for K=1, and p[1]=1p, p[0]=0, p[1]=p agrees with the output of GCFcf(N,n,p,t)
[Corrected March 30, 2020 (Thanks to Johnny Fonseca]
Hint: You still use rsolve, but now you need more initial conditions, that later on you have to solve for using
the boundary conditions on the other hand. The order of the recurrence is 2*K, the initial conditions
are x[(K1)]=1, x[1]=1, x[0]=1, x[1]=A[1], ..., x[K]=A[K], where the A[1], ..., A[k] are to be determined.
Then you use rsolve to get a monster expression featuring the indexed variables p[(K1)], ..., p[K] AND
A[1], ..., A[K]. Then you use the fact that
x(N)=1, x(N+1)=1, ..., x(N+K1)=1, get a system of equations for {A[1], ..., A[K]} solve them,
get an answer, then substitute this to the monster expression.

Write a procedure
ProbWinningGraph(A,SetOfLosingSites,SetOfWinningSites)
that inputs a matrix A of numbers (the transition matrix) given as a list of lists, such that
A[i][j] is the probability that if today are at i, then tomorrow you would be at j
(of course it has to be a stochastic matrix, each row should have nonnegative entries, and they
should addup to 1). The walk is on that graph with vertices {1, ..., N=nops(A)}. StartingPoint is location
from 1 to N, and SetOfLosingSites and SetOfWinningSites are disjoint subsets of {1,..,N}.
The rows correpsponding to members of SetOfLosingSites and SetOfWinningSites should be unit vectors
that definitely point to the same place.
For example, for the simple gambler's ruin problem with N=3 the matrix is (calling n=0 1 to accomodate
Maple's convention of indexing)
A=[[1,0,0,0],[1p,0,p,0],[0,1p,0,p],[0,0,0,1]]
SetOfLosingSites={1} SetOfLosingWinning={4}
The output should be a vector of length N that tells you the probability of winning at each starting point.
Of course the locations occupied by SetOfLosingSites is automatically 0 and
the locations occupied by SetOfWinningSites is automatically 1

Write a procedure
ExpDurationGraph(A,SetOfLosingSites,SetOfWinningSites)
that does the analogous thing for expected duration.
Of course the locations occupied by SetOfLosingSites and
by SetOfWinningSites is automatically 0

Write a procedure
PGFGraph(A,SetOfLosingSites,SetOfWinningSites,t)
that does the analogous thing for probability generating functions for the duration starting at any location.
Of course the locations occupied by SetOfLosingSites and
by SetOfWinningSites is automatically t^{0}=1

Write a procedure
LineGraph(N,p)
that inputs a positive integer N and a number of symbol p, and finds the N+1 by N+1 matrix
correponding to the simple gramber's ruin game with maximal capital N. For example
LineGraph(3,p) should be the [[1,0,0,0],[1p,0,p,0],[0,1p,0,p],[0,0,0,1]].
Check the above three procedures with the procedures in C17.txt
Added April 6, 2020: Read the students' nice solutions.
Programs discussed during the reomote class of Thurs,, April 2, 2020
C19.txt
Contains

RandDeck(n): A random starting position with cards 1, ..., 2*n, The output
is a list of length 2 with the decks of the first player

OneMove(P): One move in a 2player War where the winning card and card won goes to the BOTTOM of the winner's deck
with the won card exactly at the bottom, and the winning card on top of it. For example
OneMove([[3,1,4],[2,5,6]]); should return [[1,4,3,2],[5,6]].

OneGameV(n,K): VERBOSE VERSION.Simulates one game with 2*n cards, two players and the deck being {1, ..., 2*n}
it prints out all the intermedia steps. K is the maximum length allowed. Try:
OneGameV(4,100);

OneGame(n,K): TERSE version.Simulates one game with 2*n cards, two players and the deck being {1, ..., 2*n}
Returns the pair [winner, numer of moves]. Try:
OneGame(4,100);

ManyGames(n,K,t,M): inputs a positive integer n and a large pos. integer K, a variable t and another large positive integer M, runs
outputs the triple
[G.f for duration if Left won,G.f for duration if Right won,Number of Games that did not end in K moves]
Try:
ManyGames(3,100,t,1000);

DistD(P): The distance of a pair from the end. Try:
DistD([[1,2,3],[4,5,6]]);

WarGF(n,t): The distance generating function to the end of all initial decks in a 2player War game with deck {1, ..., 2*n}. Try:
WarGF(3,t);

[Added April 2, 2020 before class]:
WarGFplus(n,t): The distance generating function to the end of all initial decks
also returns the maximal length and the champions. Try:
WarGFplus(4,t);

[Added April 2, 2020 before class]:
PlayGame(INI):Plays a game with initial deck INI. Try:
PlayGame([[1,2],[3,4]]);

[Added April 2, 2020 before class]:
FIN(INI):Plays a game with initial deck INI. Only outputs the number of moves and the final deck. Try:
FIN([[1,2],[3,4]]);
Homework for Thurs., April 2, 2020, class (due Sunday April 5, 10:00pm)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#19
and an attachment
hw19FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

Following Yuxuan Yang's and Jingyang Deng's discoveries of periodic orbits, modify procedure
PerDecks(n)
that inputs a positive integer n and finds all the cycles in the War graph (that has (2n)! vertices).
Hopefully it would terminate with n=5, and who knows? even with n=6.
If possible, characterize these periodic orbits, and genealize them.

Fix procedure WarGF(n,t) so that it excludes members of periodic orbits, and find as many terms of
the sequence "maximal length of a game that terminates". See whether it is in the OEIS.

[Optional, if it is too much work, you can skip it. People who would pick this project would do it later]
Generalize procedures, RanDeck(n), call it GenRanDeck(n,k,r) where the deck consists of k identical copies of a deck
labeled {1, ....n,} and there are r players (so it would be a list of length r each with n*k/r cards). Then
generalize OneMove, OneGame, and ManyGames.

Get ready for the next class, by reading this
masterpiece
Added April 6, 2020: Read the students' nice solutions.
Programs discussed during the reomote class of Monday, April 6, 2020
C20.txt
See the descriptions in the file.
Reading homework for Monday April 6, 2020, class (due Thurs. April 9, 10:00pm)
Read and understand the wikipedia article on the Nash Equilibrium
Maple Homework for Monday April 6, 2020, class (due Sunday April 12, 10:00pm)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#20
and an attachment
hw20FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

A good way to describe a discrete probability distribution (especially in computer algebra)
is via the probability generating function. For example if your "Roulette" has probability
1/5 of landing on 3, Probability 1/2 of landing on 5, and probability 3/10 of landing on 20, then
the probability generating function, in the formal variable s is
1/5*s^3 + 1/2*s^5+3/10*s^(20)
Generalize the procedures , ProbWinning(L,p), ExpDuration(L,p), and PGF(L,p,t) of
C20.txt
Call them
ProbWinningG(L,s,p), ExpDurationG(L,s,p), and PGFG(L,s,p,t) ,
where now L is no longer a list of integers between 1 and min(i,Ni) but a list of probability
generating functions where L[i] is the probability generating function of the stake that you would
make if you had i dollars. In particular, for checking purposes,
ProbWinningG([s$(N1)],s,p); should give the same as ProbWinning(TiS(N),p);
and
ProbWinningG([seq(s^min(i,Ni),i=1..N1)],s,p); should give the same as ProbWinning(BoS(N),p);
Comment: The team whose project would be to explore Strategic gambling (if it exists) would
experiment with this.

Generalize GRsim1(L,p,n) , GRsim(L,p,n,M) likewise. You may use procedure Roulette(L) from
C16.txt .
Added April 13, 2020: Read the students' nice solutions.
Programs discussed during the reomote class of Thurs., April 9, 2020
C21.txt
See the descriptions in the file.
Maple Homework for Thurs. April 9, 2020, class (due Sunday April 12, 10:00pm)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#21
and an attachment
hw21FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

Write procedures BRmixedRow(G,q) and BRmixedCol(G,p) for the set of Best Responses (now possibly infinite subses of the
unit inerval) of The Row player to the mixed strategy [q,1q] of the Column Player, an
Column Player to the mixed stragegy of the Row player for a twoplayer game with two strategies available for each player

[A little bit of a challenge] Use these to write a procedure
MNE(G)
that inputs a 2 by 2 bimatrix G and outputs the set of mixed Nash Equilibria.

[Optional Challenge] Write a procedure
PNE3(G)
That inputs a threedimensional array describing a game with THREE players, where player ONE
has A strategies, Player TWO has B strategies and Player THREE has C strategies
and G is a listoflistsoflists such that G[i][j][k] is the triple [PayOffForOne,PayOffForTwo,PayOffForThree]
if they adopt strategies [i,j,k] respectively.
What is PNE3(G) if
G:=
[
[ [[3,3,3],[0,4,2]], [[4,0,1],[1,1,2]] ],
[ [[1,5,3],[2,1,3]], [[2,0,1],[1,1,2]] ]
]:
Also write a procedure
RandGameThree(A,B,C,K1,K2): A random 3player game
with A available strategies for Player ONE
and B available strategies for Player TWO
and C available strategies for Player TWO
with payoffs between K1 and K2. Then try out PNE3(G) on many such random games.

[Optional (possibly difficult) Challenge]
Write MNE(G) for an arbitrary twoperson game, or at least for the 3 by 3 case.
Added April 13, 2020: Read the students' nice solutions.
Programs discussed during the reomote class of Monday, April 13, 2020
C22.txt
See the descriptions in the file.
Maple Homework for Monday April 13, 2020, class (due Sunday April 19, 10:00pm)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#22
and an attachment
hw22FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

By using procedures RandGameG(L,K1,K2), and PNEg(G), in C22.txt ,
write a procedure
EstDist(L,K1,K2,M)
that inputs a list L describing the fact that there are n:=nops(L) players, player i (from 1 to n) has
L[i] possible strategies and the payoffs are between K1 and K2 (later make K1=0 and K2 VERY big, say 100000000000).
It also outputs a large integer M. It should output an approximate distribution for the number of
Pure Nash Equilbiria in a game with description L, with coefficeints in floating point.
Please record (we will compare notes)
EstDist(L,0,10^10,10000)
for L=[2,2],L=[3,3], L=[10,10], L=[2,2,2], L=[2,2,2,2], and L=[2,2,2,2,2]
If your computer would allow, you are welcome to replace 10000 by 100000.

Write a procedure
Centepede(L)
That inputs a list of odd length, call it 2n+1, describing the payoffs for Player I and Player II
(in that order) in an ndouble round Centepeded game, where the last two entries desribe the terminal
payoffs for player II, and (2*i1)th entry is the payoff for Player I if he quits at the ith doubleround
and the 2ith entry is the payoff for Player II.
It outputs the payoffs under socalled "rational play"
using backwards induction. For exaple, in the game described in Appendix 2 of this
interesting opinion
Centepede([[2,1],[1,2],[4,1],[6,3],[3,4]]);
should output [2,1].

Using Centepede(L), write a program
Contract(L) ,
that decides whether it is good for the two players sign a contract
not to use one of the terminal options of player II, because it would benefit both of them.
It should output the modified (improved) payoff pairs. For example
Contract([[2,1],[1,2],[4,1],[6,3],[3,4]]);
should return
true,[6,3]
Added April 20, 2020: Read the students' nice solutions.
Programs discussed during the reomote class of Thurs., April 16, 2020
Maple Homework for Thurs. April 16, 2020, class (due Sunday April 19, 10:00pm)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#23
and an attachment
hw23FirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

Modify Alison Bu's beautiful procedure E(m,p) in
C23.txt
to write a procedure
PGF(m,p,t)
that inputs positive integers m and p and a variable name t, and outputs the probability generating
function for the random variable "duration of game"

Using PGF(m,p,t) and what we had previously, write a procedure
Moms(m,p,K)
that inputs m and p as above and outputs a list of length K whose first entry is the expectation,
second entry is the standarddeviation, and the ith entry is the ith scaled moment about the mean.
What are
Moms(100,110,6); Moms(200,220,6); Moms(300,660,6);
Do the third through sixth's entry agree with each other?

By playing around with Alison's procedure, E(m,p), I observerd that for fixed k
E(m,m+k) for large m grows like a constant C(k)*sqrt(m). Write a procedure
EstimateCk(k), that estimates this constant (if it exists)

[Optional challenge, 25 dollars donation to the OEIS in honor of the first provers. I do not know the answer]
Prove the above , and if possible find an explicit expression for C(k) in terms of k.

[Optional challenge, 25 dollars donation to the OEIS in honor of the first provers. I do not know the answer,
or whether there is an answer]
Conjecture and prove a closedform formula for E(m,p).
Added April 20, 2020: Read the students' nice solutions.
Programs discussed during the reomote class of Mon., April 20, 2020
Reading/Listening Homework for Monday April 20, 2020, class (due Thurs. April 23, 11:00am)
This homework assignment is in memory of my great hero John Horton Conway (19382020).
Optional Maple Homework for Monday April 20, 2020, class (due Thurs. April 23, 11:00am)
Please email ShaloshBEkhad at gmail dot com an email with
Subject: hw#24a
and an attachment
hw24aFirstNameLastName.txt
and indicate whether it is OK to post. If you do not say "Please do not post", it will be posted.

By reading Dr. Z.'s brief introduction to
Conway's Surreal numbers (or wikipedia, or Knuth's book, or wherever you can find it),
write a Maple "Surreal Numbers Calculators". Note that evetyhing is RECURSIVE, so it may be a good idea to
have "option rememeber".
It should contain at least the following procedures

[Super optional] Program the "Look and Say Sequence" in Latin Numerals (see the end of the above mentioned
video). Determine its rate of growth, empirically
finding the "Latin Numeral Conway constant".
Homework for Monday April 20, 2020, class (due Sunday. April 26, 10:00pm)

Implement an analog of the epidemic model
Flahaut et. al.'s , paper
for two cities (and still four types), getting a system of 8 differential equations.
It should have as, input a time step, h (supposed to be small), and the derivatives f'(t) should be
approximated by (f(t+h)f(t))/h. See whether you can similar graphs for the growth of the epidemics.

[Optional] Implement the above paper
for an arbitrary number,n, of cities (not just 2 or 9), getting a system of 4n differential equations.
If course the "transition matrix" c_{ij} should be part of the input. To test it
write a program to generate random n by n such transition matrices.

[Super Opptional, no deadline]
Continue the preliminary work in
C24.txt
implementing Hurd's paper.
Added April 27, 2020: Read the students' nice solutions.
Programs discussed during the reomote class of Thurs, April 23, 2020

C25.txt : A Surreal number calculator.
Optional Homework for Thurs. April 23, 2020, class (due Sunday. April 26, 10:00pm)
Write a procedure MUL(X,Y) that multiplies two surreal numbers X and Y
Added April 27, 2020: Read Charles Kenny's nice solution
to the optional problem to multiply two surreal numbers, now also part of
C25.txt .
Programs discussed during the reomote class of Monday April 27, 2020

C26.txt : Exact answers (but very slow, so far for,
for more than two players) for the distributions of the
number of Pure Nash Equilibria in nplayer games with an arbitrary (but finite) number of strategy choices for each player.
Contains procedures:

RandGameGm(L): A random game with n:=nops(L) players where each of the players has distinct payoffs equal to
{1, ..., L[1]*...*L[n]}

IsComp(S): inputs a set of strategies, S, in an nplayer game outputs true iff they are compatible (see definition below)
Definition: A set S of r strategy lists in an nplayer game
S={[m1[1], ..., m1[n]], [m2[1], ..., m2[n]] , [m3[1], ..., m3[n]], ..., [mr[1], ..., mr[n]], ..., }
is compatible if for each i from 1 to n the set of vectors of length n1
S_{i}={
[m1[1], ..., m1[i1],m1[i+1], ..., m1[n] ],
[m2[1], ..., m2[i1],m2[i+1], ..., m2[n] ],
...
[mr[1], ..., mr[i1],mr[i+1], ..., mr[n] ]
}
has the same cardinality as S, i.e. r. This is because for each choicevector of the other n1 players
the best response of player i should be unique (since all its payoffs are distinct).

WtSeq(L): The weight sequence if the list of available strategies is L, it is a list of length mul(L[i],i=1..nops(L))
(but the larger entries are 0) such that the ith entry is the sum of the weights of all
subsets of the set of Strategy vectors, where the weight of a set of strategy choices is the probability
that they are all Nash Equilibria. It is 0 if the set is not compatible (see above) but otherwise
it is 1/mul(L[i],i=1..nops(L)) to the power its cardinality.
Try:
WtSeq([2,2]);

ZLs(L,s): The prob. that there are exactly s Pure Nash Equilibria in a random nops(L)player game where Players i has L[i] strategy choices.
ZLs([2,2],0);

ZLt(L,t): The polynomial in t whose coefficient of t^s is the prob. of having exactly s pure Nash Equilibria in
a random nops(L)player game where Player i has L[i] strategy choices. Warning, VERY slow (but exact)
and for each player all its payoffs are distinct.
It is the exact value of what the following line approximates
add(t^nops(PNEg(RandGameGm(L))),i=1..10000)/10000.;
Try:
ZLt([2,2,2],t);

Zmns(m,n,s): The prob. that there are exactly s Pure Nash Equilibria in a random 2player game where Players 1 and 2 having m and n strategies,
using the exact formula that employs the Principle of Inclusion Exclusion and the fact that
the number of compatible subsets of {[i,j]; 1 ≤ i ≤ n, 1 ≤ j ≤ m} with s members (s between 0 and min(m,n)) is
exactly binomial(n,r)^2*r! [Writing the members in lexicographic order, there are binomial(n,r) choices for the
top row (regarding Player 1), and n(n1)*...*(nr+1)=binomial(n,r)*r! choices for the second row (regarding player 2).
It is an efficient way to compute ZLs([m,n],s) (see above) for the special case of two players.

Zmnt(m,n,t): The polynomial of degree min(m,n) in t whose coefficient of t^s is the prob. of having exactly s pure Nash Equilibria in
a random 2player game where Players 1 and 2 having m and n strategies
It is an efficient way to compute ZLt([m,n],t) (see above) for the special case of two players.
Optional Challenging Homework for Monday April 27, 2020, class (no due date)
If and when you find an answer email me at DoronZeil at gmail dot com (NOT ShaloshBEkhad at gmail dot com)

[20 dollars donation to the OEIS in honor of the prover]
Accodring to the above procedure, Zmns(m,n,s), for the special case m=n,
the EXACT probability that in a 2player game where
each player has n strategy choices is given by
sum((1)^(rs)*binomial(r,s)*binomial(n,r)^2*r!/n^(2*r),r=s..n)
Prove that (as is apparently wellknown, but possibly via another method) fact, that as n goes to infinity, and s fixed,
this converges to 1/(e*s!). [If you can do the special case s=0, I will pay 10 dollars].
[Added April 27, 2020, 4:35pm, after talking to RDB: It may not be as hard as I thoght, for fixed r and large n
and then as n goes to infinity
binomial(n,r)*r! is roughly n^r,
one gets 11+1/2!1/3!+... =1/e, but to make it fully rigorous may take some effort.

[50 dollars donation to the OEIS in honor of the prover] Find an EFFICIENT way to compute
ZLs(L,s) ,
at least for the threeplayer case, using (if possible) the Principle of InclusionExclusion.
Comment: For the twoplayer case the problem of finding the number of compatible sets of a given cardianlity,
namely, WtSeq([m,n]) was almost trivial (see above). Finding WtSeq([m1,m2,m3]) seems much harder
(unless I missed something).
Programs discussed during the reomote class of Thurs. April 30, 2020
Today we discussed this
great paper.
There is no homework, except that hopefully the "Gambler's Ruin" team will extend it. See the notes in the
project page.
Optional Challenging Homework for Thurs. April 30, 2020, due Monday, May 4 during class
Prizes for the best entries

In yesterday's New Your Times (April 29, 2020) there was the following puzzle
Find three letters such that if you prefix it to "eli" and append them in reverse order, you would
a common phrase consisting of two words (without a space between them).
Answer: [t,i,m]. Indeed the reverse of [t,i,m] is [m,i,t] and you get
time limit
Using procedure GREP(N,S) in
C28.txt
write a Maple program that
(i) Solves such puzzles
(ii) Makes them up
Added May 4, 2020: See the nice Maple Code solving this problem by
Chun "Larry" Lau
and Yukun Yao

In today's New Your Times (April 30, 2020) there was the following (classical CRYPTOGRAM) puzzle
Find a permutation of the alphabet (26! choices if done by complete brute force) such
that the following sentence makes sense [capitalization ignored, using all lowercase].
In other words solve the CRYPTOGRAM
wh wm r pjwzc fk hdx xlsewmd erlsjrsx hdrh yfhdxz rly
yrex mdrzx r kwzmh exhhxz rm vf krhdz rlv kxyrex
Using procedure GREP(N,S) in
C28.txt
write a Maple proedure that
(i) Solves such puzzles. It would output the permutation of [a, ..., z] and the
deciphered sentence
(ii) Makes them up (starting with an English sentence
Notes:
1.input should be formatted to conform to the convention of GREP(N,S), i.e.
a list of lists, so the above puzzle would become
[[w,h], [w,m], [r], [p,j,w,z,c], [f,k], ..., [k, x, y,r, e, x]]
2. DO NOT use lowercase letters for programming, and indexing, these would conflict with the list called ENGLISH in
C28.txt. Use CAPIAL LETTERS for Maple programming: e.g. "for A from 1 to nops(ENGLISH) do"
Programs discussed during the reomote class of Monday, May 4, 2020
Class Projects

The War Project. Team leader: Robert DoughertyBliss. Team members: Charles Kenny, Alison Bu
.pdf

Generalized Gambler's Ruin (Leader: Yukun Yao; Team members: Victoria Chayes, Chun Lau, Yuxuan Yang)
Here are the current version of the
paper and
the
Maple program.

The Nash Equilibria project (Leader: Alan Chenoff; Team members: Khizar Anjun, Jingyang Deng, Zidong Zhang)
Here is the current version of this project, Maple package
EM20ProjNash.txt.

Strict War: (Team leader and only member, Johnny Fonsenca)
Added May 20, 2020: Read Johnny Fonseca's
interesting article
about Strict War.
Dr. Z.'s teaching page