###################################################################### ## MVNM.txt Save this file as MVNM.txt to use it, # # stay in the # ## same directory, get into Maple (by typing: maple ) # ## and then type: read `MVNM.txt` # ## Then follow the instructions given there # ## # ## Written by Doron Zeilberger, Rutgers University , # ## DoronZeil at gmail dot com # ###################################################################### print(` Written: Jan./Feb. 2022: tested for Maple 2020 `): print(): print(`This is MVNM.txt, A Maple package`): print(`accompanying Shalsoh B. Ekhad and Doron Zeilberger's article: `): print(` A Linear Time, and Constant Space, Algorithm to Compute the Mixed Moments of the Mutlivariate Normal Distribution `): print(): print(`The most current version is available on WWW at:`): print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/MVNM.txt .`): print(`Please report all bugs to: DoronZeil at gmail dot com .`): print(): print(`For general help, and a list of the MAIN functions,`): print(` type "ezra();". For specific help type "ezra(procedure_name);" `): print(`For a list of the supporting functions type: ezra1();`): print(): ezra1:=proc() if args=NULL then print(`The SUPPORTING or Deprated procedures are`): print(` MOM2d, MOM3d`): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(` MVNM.txt : A Maple package for efficient computation of mixed moments of Bivariate and Trivariate Normal Distribution`): print(`The MAIN procedures are: MOM2, MOM3, MOMd`): print(``): elif nargs=1 and args[1]=MOM2 then print(`MOM2(c,m0) Computes Int(int(x^m0[1]*y^m0[2]*F(x,y,c),x=-infinity..infinity), y=-infinity..infinity) using the REURRENCES. c can be symboic or numeric`): print(`Try:`): print(`MOM2(c,[5,5]);`): elif nargs=1 and args[1]=MOM2d then print(`MOM2d(c,m0) Computes Int(int(x^m0[1]*y^m0[2]*F(x,y,c),x=-infinity..infinity), y=-infinity..infinity) DIRECTLY using the moment generationg function`): print(`Try:`): print(`MOM2d(1/2,[5,5]);`): elif nargs=1 and args[1]=MOM3 then print(`MOM3(c12,c13,c23,m0): The (m0[1],m0[2],m0[3]) mixed moment of the trivariate normal distribution with Covariance matrix.`): print(`[[1,c12,c13],[c12,1,c23],[c13,c23,1]].`): print(` Done via the RECURRENCES. c12,c13,c23, can be numeric or symbolic, Try:`): print(`MOM3(c12,c13,c23,[10,10,10]);`): elif nargs=1 and args[1]=MOM3d then print(`MOM3d(c12,c13,c23,m0): The (m0[1],m0[2],m0[3]) mixed moment of the trivariate normal distribution with Covariance matrix`): print(`[[1,c12,c13],[c12,1,c23],[c13,c23,1]]. Try:`): print(`, using the moment generating function. Try`): print(`MOM3d(c12,c13,c23,[10,10,10]);`): elif nargs=1 and args[1]=MOMd then print(`MOMd(C,m): Given a covariance matrix, say of dimension n by n, and a list m of length n, outputs`): print(`the mixed -m moments of the multivariate normal distribution (with 0 mean and unit variancs) with that covariance matrix. `): print(`C can be numeric or symbolic. Try:`): print(`MOMd([[1,c],[c,1]],[6,8]);`): else print(`There is no such thing as`, args): fi: end: with(linalg): MOM2d:=proc(c,m0) local gu,t1,t2,z: if m0[1]+m0[2] mod 2 =1 then RETURN(0): fi: gu:=exp((1/2*t1^2+c*t1*t2+1/2*t2^2)*z): gu:=taylor(gu,z=0,(m0[1]+m0[2])/2+1): gu:=expand(coeff(gu,z,(m0[1]+m0[2])/2)): m0[1]!*m0[2]!*coeff(coeff(gu,t1,m0[1]),t2,m0[2]): end: MOM3d:=proc(c12,c13,c23,m0) local gu,t1,t2,z: if m0[1]+m0[2]+m0[3] mod 2 =1 then RETURN(0): fi: gu:=exp((1/2*t1^2+1/2*t2^2+1/2*t3^2+c12*t1*t2+c13*t1*t3+c23*t2*t3)*z): gu:=taylor(gu,z=0,(m0[1]+m0[2]+m0[3])/2+1): gu:=expand(coeff(gu,z,(m0[1]+m0[2]+m0[3])/2)): m0[1]!*m0[2]!*m0[3]!*coeff(coeff(coeff(gu,t1,m0[1]),t2,m0[2]),t3,m0[3]): end: EvalScheme1:=proc(S,m,M,m0) local ope,INI,L,a,kha,j: ope:=S[1][1]: INI:=S[2]: L:=INI: if m0<=nops(INI)-1 then RETURN(INI[m0+1]): fi: for a from nops(INI) to m0 do kha:=normal(add(subs(m=a,coeff(ope,M,-j))*L[-j],j=1..-ldegree(ope,M))): L:=[op(2..nops(L),L),kha]: od: L[-1]: end: #EvalScheme(S,m,M,m0): inputs a scheme S (let n:=nops(S)) in n dimension, phrased in terms of m and M (i.e. the discrete variables are m[1],...,m[n] and the corresponding shift operators are #M[1],..., M[n], and a list of non-negative integers m0. Evaluats it at m0. Try #EvalScheme([[2/M[1],3/M[2],4/M[3]],[[[1]]],m,M,[2,2,2]); #EvalScheme([[m[1]/M[1]],[1]],m,M,[3]); EvalScheme:=proc(S,m,M,m0) local n,ope1,GU,S1,A1,INI1,a,i1,m01,i: n:=nops(S[1]): if nops(m0)<>n then RETURN(FAIL): fi: if n=1 then RETURN(EvalScheme1(S,m[1],M[1],m0[1])): fi: ope1:=S[1][1]: A1:=-ldegree(ope1,M[1]): if m0[1]