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: