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

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 computer-algebra 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 golden-oldie 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

Programs done on Thurs., Jan. 23, 2020

C1.txt , Contains procedures

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.

  1. [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 .

  2. [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.

  3. [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

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.

  1. [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 30-60 of Frank Garvan's awesome Maple booklet ( part 1, part 2)
    Only record a small sample in hw2YourName.txt .

  2. Read and understand Chapter 2, sections 2.1-2.4 (pp. 17-22) of Hoff's book.

  3. 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).

  4. 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 simple-minded 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 .

  5. [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 N-k Tails) is

    N!/(k!(N-k)!)pk(1-p)N-k

Added Feb. 2, 2020: Read the students' nice solutions

Programs done on Thurs., Jan. 30, 2020

C3.txt , Contains procedures

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.

  1. [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 60-90 of Frank Garvan's awesome Maple booklet ( part 1, part 2)
    Only record a small sample in hw3YourName.txt .

  2. Read and understand this this masterpiece

  3. 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.

  4. We will soon learn that the PGF of the binomial distribution is (p*t+(1-p))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 3-rd through 12-th 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

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.

  1. [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 .

  2. 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*(1-p)^(n-k),k=0..trunc(n*p+x*sqrt(n*p*(1-p)))):
    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?

  3. 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*(1-a/n)^(n-k),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.

  4. [Human exercise]: Prove it humanly.

  5. [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 On-Line Encyclopedia of Integer Sequences   .

Optional Homework for Thurs., Feb. 6, 2020, class (due indefnitely)

  1. [$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

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.

  1. Write a procedure:

    CheckHT(n,p,x,N), that first calls HT(n,p,x) to find out what is the "cut-off" below which we have to reject the "null hypothesis" that the coin has Pr(Heads)=p with (one-sided) 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:

  2. 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*(1-p)^k

  3. Read section 3.1 of Hoff's book

  4. 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/2-z,1/2+z], that contains a portion x of the distribution, i.e. such that

    int(f,t=1/2-z..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)

  5. [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)

  6. [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

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.


  1. [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 Dougherty-Bliss 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.

  1. Read Chapter 4 of Hoff's book

  2. 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.

  3. 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 Dougherty-Bliss . Then we explored some of the built-in 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.

  1. 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 closed-form, or its inverse function of the CDF is not closed form it should return FAIL.

  2. 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)  

  3. 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((x-a1)^2/(2*b1)*exp((y-a2)^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.

  4. 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.

  5. [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.

  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

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.

  1. [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?

  2. Wtite a continuous analog of EAv(X,f,x), call it CAv(P,t,f,a,b) that inputs a NON-normalized 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

  3. [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 "data-collecting" period of length M2, to esimate the exact value of CAv(P,t,f,a,b).

  4. 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:

Reading homework for Monday Feb. 24, 2020, class (due Thurs. Feb. 27, 2020)

  1. 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.

  1. [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?

  2. [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]

  3. 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 bulid-up, 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:

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.

  1. Implement the NON-SAMPLING 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.

  2. Extend the programs in C11.txt to the situation where there are two cross-overs, 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.

  3. [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 "floating-point", so you need to convert it into an integer by doing trunc().]

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

  5. [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:

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.

  1. 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 j-th entry is the j-th 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 j-th entry is the j-th 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 j-th entry is the j-th 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((x-mu)^2*f(x),x), and more generally, the k-th moment about the mean is Int((x-mu)^k*f(x),x). Analogously for the discrete case]

  2. 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)

  3. Do the Poisson analog of what was done in C12.txt , i.e. Gibbs Sampling for the cross-over 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 ,

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.

  1. 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

  2. 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.

  3. [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 lists-of-lists-of-lists-of-lists, ... nested k times. for example one member of ArrayS({0,1,2},[2,2,2]) is [ [[0,2],[1,2] ], [ [1,1],[2,2] ] ].

  4. 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 ,

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.

  1. 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 "wrapped-around" M by N by K box.

  2. Wtite a procedure

    NewA3d(A,c)

    that inputs a three-dimensional 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

  3. 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 ,

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.

  1. 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)

  2. Now try:

    NickAv(40,40, 1.3 , 1.0, 1000,10000, E, E,M);

    ten times. Are the outputs closer to each other?

  3. According to Julia Yeomans's book starting with either the all 1-matrix 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.

  4. Write a procdeure

    Shrink(A)

    that starts with a 3k by 3k matrix of {1,-1} and outputs a 3k-1 by 3k-1 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

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.

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

  2. 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.

  3. 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 non-positive amount or an amount >=N then the game is over

    [Note that GenGRsim1(N,[[p,1],[1-p,-1]],n) is equivalent to GR1sim(N,p,n)]

  4. Run

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

    For

    Do you notice considerable difference between the first and the second outcome? Between the second and the third outcome? Why?

  5. 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

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.

  1. 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 third-moment about the mean), the kurtosis (scaled fourth-moment about the mean)?

  2. [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 add-up to 1, but Maple does not have to be notified). The game ends as soon as the capical is non-positive or strictly larger than N-1. [Of course "winning" a negative amount is losing].

    Check that your output for K=1, and p[-1]=1-p, 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[-(K-1)]=1 and the boundary conditions are x[N]=1, x[N+1]=1, ..., x[N+K-1]=1. (Since t0=1)

Added March 29, 2020: Read the students' nice solutions.

Programs discussed during the reomote class of Mon., March 30, 2020

C18.txt

Contains

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.

  1. [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 add-up to 1, but Maple does not have to be notified). The game ends as soon as the capical is non-positive or strictly larger than N-1. [Of course "winning" a negative amount is losing].

    Check that your output for K=1, and p[-1]=1-p, 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[-(K-1)]=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[-(K-1)], ..., p[K] AND A[1], ..., A[K]. Then you use the fact that x(N)=1, x(N+1)=1, ..., x(N+K-1)=1, get a system of equations for {A[1], ..., A[K]} solve them, get an answer, then substitute this to the monster expression.

  2. 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 non-negative entries, and they should add-up 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],[1-p,0,p,0],[0,1-p,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

  3. 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

  4. 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 t0=1

  5. 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],[1-p,0,p,0],[0,1-p,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

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.

  1. 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.

  2. 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.

  3. [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.

  4. 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.

  1. 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,N-i) 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$(N-1)],s,p); should give the same as ProbWinning(TiS(N),p);

    and

    ProbWinningG([seq(s^min(i,N-i),i=1..N-1)],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.

  2. 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.

  1. 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,1-q] of the Column Player, an Column Player to the mixed stragegy of the Row player for a two-player game with two strategies available for each player

  2. [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.

  3. [Optional Challenge] Write a procedure

    PNE3(G)

    That inputs a three-dimensional 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 list-of-lists-of-lists 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 3-player 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.

  4. [Optional (possibly difficult) Challenge] Write MNE(G) for an arbitrary two-person 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.

  1. 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.

  2. 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 n-double round Centepeded game, where the last two entries desribe the terminal payoffs for player II, and (2*i-1)-th entry is the payoff for Player I if he quits at the i-th double-round and the 2i-th entry is the payoff for Player II. It outputs the payoffs under so-called "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].

  3. 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.

  1. 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"

  2. 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 standard-deviation, and the i-th entry is the i-th 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?

  3. 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)

  4. [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.

  5. [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 closed-form 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 (1938-2020).

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.

Homework for Monday April 20, 2020, class (due Sunday. April 26, 10:00pm)

  1. 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.

  2. [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" cij should be part of the input. To test it write a program to generate random n by n such transition matrices.

  3. [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

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

Contains procedures:

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)
  1. [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 2-player game where each player has n strategy choices is given by

    sum((-1)^(r-s)*binomial(r,s)*binomial(n,r)^2*r!/n^(2*r),r=s..n)

    Prove that (as is apparently well-known, 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 1-1+1/2!-1/3!+... =1/e, but to make it fully rigorous may take some effort.

  2. [50 dollars donation to the OEIS in honor of the prover] Find an EFFICIENT way to compute

    ZLs(L,s)   ,

    at least for the three-player case, using (if possible) the Principle of Inclusion-Exclusion.

    Comment: For the two-player 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

  1. 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

  2. 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 lower-case]. 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 lower-case 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

Dr. Z.'s teaching page