#C9.txt: Feb. 20, 2020
Help:=proc(): print(` RRV(K1,K2a,k2b,N), EAv(X,f,x) , DAv(X,f,x,M)`): 
print(`NickAv(X,f,x,M1,M2),  NickAvWrong(X,f,x,M1,M2)`):
end:

with(Statistics):

#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 ], ...]
RRV:=proc(K1,K2a,K2b,N) local P,ra1,ra2,i:
ra1:=rand(1..K1):
P:=[seq( ra1(), i=1..N)]:
P:=P/add(P[i],i=1..nops(P)):
ra2:=rand(K2a..K2b):
[seq([ra2(),P[i]],i=1..N)]:
end:

#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)]
EAv:=proc(X,f,x) local i:
add(subs(x=X[i][1],f)*X[i][2],i=1..nops(X)):
end:

#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")
DAv:=proc(X,f,x,M) local N,ra, T, B,i,pe:
N:=nops(X):
ra:=rand(1..N):
T:=0:
B:=0:

for i from 1 to M do
 pe:=ra():
 T:=T+subs(x=X[pe][1],f)*X[pe][2]:
 B:=B+X[pe][2]:
od:

evalf(T/B):
end:




#NickAvWrong(X,f,x,M1,M2): Original wrong version of NickAv(X,f,x,M1,M2) with the inequality r<rat reversed
#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)
NickAvWrong:=proc(X,f,x,M1,M2) local N,pe1,pe2,i,rat,U,ra,r,T:
U:=RandomVariable(Uniform(0,1)):
N:=nops(X):
ra:=rand(1..N):

pe1:=ra():

for i from 1 to M1 do
 pe2:=ra():
 rat:=X[pe2][2]/X[pe1][2]:
 if rat>=1 then
  pe1:=pe2:
  else
  r:=Sample(U,1)[1]:
  if rat<r then
   pe1:=pe2:
  fi:
 fi:
od:
  
T:=0:

for i from 1 to M2 do
 pe2:=ra():
 rat:=X[pe2][2]/X[pe1][2]:
 if rat>=1 then
  pe1:=pe2:
  else
  r:=Sample(U,1)[1]:
#The line below is erroneous. It should be r<rat
  if r>rat then
   pe1:=pe2:
  fi:
 fi:

 T:=T+subs(x=X[pe1][1],f):
od:

 evalf(T/M2):
end:



#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)
NickAv:=proc(X,f,x,M1,M2) local N,pe1,pe2,i,rat,U,ra,r,T:
U:=RandomVariable(Uniform(0,1)):
N:=nops(X):
ra:=rand(1..N):

pe1:=ra():

for i from 1 to M1 do
 pe2:=ra():
 rat:=X[pe2][2]/X[pe1][2]:
 if rat>=1 then
  pe1:=pe2:
  else
  r:=Sample(U,1)[1]:
  if rat<r then
   pe1:=pe2:
  fi:
 fi:
od:
  
T:=0:

for i from 1 to M2 do
 pe2:=ra():
 rat:=X[pe2][2]/X[pe1][2]:
 if rat>=1 then
  pe1:=pe2:
  else
  r:=Sample(U,1)[1]:
#correction after class, thanks to Jingyang Deng
  if r<rat then
   pe1:=pe2:
  fi:
 fi:

 T:=T+subs(x=X[pe1][1],f):
od:

 evalf(T/M2):
end: