#OK to post homework
#Jingyang Deng, Jan. 27, 2020, Assignment 2

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

Help:=proc():print(`Nor(L), LD(L), RSg(N,L), Plt(L), RSgM(N,L), Plt_by_Victoria(L)`):end:

#Normalize the list, L
Nor:=proc(L):
if min(L)<0 then
  RETURN(FAIL):
elif add(L)=0 then
  RETURN(FAIL):
else
  L/add(L):
fi:
end:

#generate a loaded die with the prob. distrib. of Nor(L)
LD:=proc(L) local L1,i,j,r,N:
L1:=Nor(L):
L1:=lcm(op(denom(L1)))*L1:
N:=add(L1):
r:=rand(1..N)():
#print(r):
for i from 1 to nops(L1) while r>sum(L1[j],j=1..i) do
od:
i:
end:

#generate a random sequence that records the results of rolling the loaded die above
RSg:=proc(N,L) local i:
[seq(LD(L), i=1..N)]:
end:

###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? 

#use Plt(L) to run the line: Statistics[Histogram]([seq(add(RSg(3000,L)),i=1..3000)]); 
Plt:=proc(L) local i:
Statistics[Histogram]([seq(add(RSg(3000,L)),i=1..3000)]):
end:

#the shape of four results are familiar. In fact, this can be explained by Central Limit Theorem.
#the random varible, add(RSg(3000,L)), is approximately subject to normal distribution.
###It took me very long time (about 3.6 hours) to plot those four histograms.

#the probability of an "atomic" event (the outcome of one experiment in {H,T}^n) 
#is p^(#heads)*(1-p)^(#tails), and there are n choose k, i.e., n!/(k!*(n-k)!)
#"atomic" events contained in the "compound event" of interest. 
#Hence, the prob.=n!/(k!*(n-k)!))*p^k*(1-p)^(n-k).

#[Added Jan. 28, 2020]
#Code by Victoria Chayes
with(LinearAlgebra):
with(Statistics): 
RSgM:=proc(N,L) local X:
X:=RandomVariable(ProbabilityTable(L)):
convert(Sample(X,N),list):
end:

#use Plt_by_Victoria(L) to run the line: Statistics[Histogram]([seq(add(RSgM(3000,L)),i=1..3000)]); 
Plt_by_Victoria:=proc(L) local i:
Statistics[Histogram]([seq(add(RSgM(3000,L)),i=1..3000)]):
end: