Help:=proc(): print(`Jump(i,L), MH(i,L,K) , RealAve(L,i,ex1),MHAve(L,i,ex1,B1,E1) `): end:


#Inputs a list L of size n, say, of relative importance
#(prob.) performs one step in the Metropolis-Hastings algorithm
#(it picks any other index,j, uniformly, if it is more important
#it definitely goes there, if it is less, then it goes there
#with prob. L[j]/L[i]
Jump:=proc(i,L) local j, ra, rat,x:
ra:=rand(1..nops(L)):
j:=ra():

rat:=L[j]/L[i]:

if rat>=1 then
 RETURN(j):
fi:

x:=stats[random,uniform[0,1]](1);

if x<=rat then
 RETURN(j):
else
 RETURN(i):
fi:

end:


#MH(K): the Markov Chain of length K starting  at i
MH:=proc(i,L,K) local M,c,j:
j:=i:
M:=[j]:

for c from 1 to K do
 j:=Jump(j,L):
  M:=[op(M),j]:
od:
M:

end:


#RealAve(L,i,ex1): the weighted 
#average of ex1(i) (ex1 is a symbolic function of index i)
#according to the weight Pr(i)=L[i]/convert(L,`+`)
#
RealAve:=proc(L,i,ex1) local j:

add(L[j]* subs(i=j , ex1)  , j=1..nops(L))/convert(L,`+`) :

end:


#MHAve(L,i,ex1): an appx. to the weighted average, according to L
#of ex1(j) (where ex1 is an expression in the symbol i)
#you run the MH algorithm B1 times for fun, and start
#averaging until E1
MHAve:=proc(L,i,ex1,B1,E1) local T,j,c:

j:=rand(1..nops(L))();

for c from 1 to B1 do
 j:=Jump(j,L):
od:

T:=0:

for c from B1+1 to E1 do
 j:=Jump(j,L):
 T:=T+evalf(subs(i=j,ex1)):
od:

T/(E1-B1):

end: