npde[ker] := (y,h,x,n) -> evalf( sum(exp(-((y-x[i])/h)^2/2),i=1..n)/(sqrt(2*Pi)*n*h) ); npde[gauss_kernel]:=proc (y,h,x,n) local c,s,j; c := sqrt(2.*Pi)*h*n; s := 0.; for j from 1 to n do s := s + exp(-((y-x[j])/h)^2/2); od; RETURN( s/c ); end; # npde[gauss_kernel] npde[vg2] := proc(y,p,m,s) # value of pdf p*N(m[1],s[1]^2) + (1-p)*N(m[2],s[2]^2) local a,b; a := p*exp(-((y-m[1])/s[1])^2/2)/s[1]; b := (1-p)*exp(-((y-m[2])/s[2])^2/2)/s[2]; RETURN( (a+b)/sqrt(2*Pi) ); end; npde[show] := proc(h) local y,p1,p2,p3; p1 := plot(vg2(y,p,m,s),y=-6..6): p2 := plot(ker(y,h,x,n),y=-6..6,color=blue): p3 := plot(ker(y,h_hat_approx,x,n),y=-6..6,color=green): plots[display]({p1,p2,p3}); end; npde[show_bw] := proc(h) # show in black and white local y,p1,p2,p3; # true is solid line p1 := plot(vg2(y,p,m,s),y=-6..6,linestyle=0): p2 := plot(ker(y,h,x,n),y=-6..6,color=blue,linestyle=2): p3 := plot(ker(y,h_hat_approx,x,n),y=-6..6,color=green,linestyle=3): plots[display]({p1,p2,p3}); end; npde[show2] := proc(h1,h2) local y,p1,p2,p3; p1 := plot(ker(y,h1,x,n),y=-5..5,color=blue): p2 := plot(ker(y,h2,x,n),y=-5..5,color=green): p3 := plot(exp(-y^2/2)/sqrt(2*Pi),y=-5..5): plots[display]({p1,p2,p3}); end; npde[g2]:=proc (n,p,m,s) # generates n obs. from pdf p*N(m[1],s[1]^2) + (1-p)*N(m[2],s[2]^2) local m1,m2,s1,s2,u,x; m1 := m[1]; m2 := m[2]; s1 := s[1]; s2 := s[2]; for i to n do u := stats[random,uniform[0,1]](); if u < p then x[i] := stats[random,normald[m1,s1]](); else x[i] := stats[random,normald[m2,s2]](); fi; od; RETURN(convert(x,list)); end; # g2 npde[one_iter] := proc(an,am2,x,i,n) local S,T,U,T2,new_i_k,k,lan,lam2; lan := an; lam2 := am2; for k to n do U := stats[random,uniform[0,1]](); new_i_k := floor( 1 + n*U ); while (new_i_k = k or new_i_k = i[k]) do U := stats[random,uniform[0,1]](); new_i_k := floor( 1 + n*U ); od; # print('new_i_k'=new_i_k,'k'=k, 'i[k]'=i[k]); T := lam2 + (x[new_i_k] - x[i[k]])*(x[new_i_k]+x[i[k]] - 2*x[k]); T2 := T^(-n/2); S := T2/lan; if S > 1 then lan := T2; lam2 := T; i[k] := new_i_k; else U := stats[random,uniform[0,1]](); if (U < S) then lan := T2; lam2 := T; i[k] := new_i_k; fi; fi; od; RETURN(lan,lam2,eval(i)); end; npde[EhA] :=proc (n_iter, n_burn, n_inter, x, n) local a,an,am2,i,estimate,k,L,ave,j,v; # start from i=[2,3,...,1] i := linalg[vector]([seq((j mod n) + 1, j=1..n)]); j := 'j'; am2 := sum( (x[j]-x[i[j]])^2,j=1..n ); an := am2^(-n/2); a := evalf( GAMMA((n-1)/2)/GAMMA(n/2)/sqrt(2) ); # burn the first n_burn for k to n_burn do L := one_iter(an,am2,x,i,n); an := L[1]; am2 := L[2]; od; # compute the average sequentially ave := 0; for k to n_iter do # skip n_inter iterations for j to n_inter do L := one_iter(an,am2,x,i,n); an := L[1]; am2 := L[2]; od; ave := ave+sqrt(am2); v[k] := a*ave/k; printf(`v[%d]=%f `,k,v[k]); od; RETURN( convert(v,array) ); end; # EhA npde[Eh] :=proc (n_iter, n_burn, x, n) local a,an,am2,i,estimate,k,L,ave,j; i := linalg[vector]([seq((j mod n) + 1, j=1..n)]); print('ok0'); j := 'j'; am2 := sum( (x[j]-x[i[j]])^2,j=1..n ); an := am2^(-n/2); print('ok1'); for k to n_burn do L := one_iter(an,am2,x,i,n); an := L[1]; am2 := L[2]; od; print('ok2'); ave := 0; for k to n_iter do L := one_iter(an,am2,x,i,n); an := L[1]; am2 := L[2]; ave := ave+sqrt(am2); od; a := evalf( GAMMA((n-1)/2)/GAMMA(n/2)/sqrt(2) ); RETURN( a*ave/n_iter ); end; # Eh npde[exact_Eh] :=proc (x,n) local i,a,j,k,N,am2,s1,s2; for j to n do i[j] := 1; od; i[1] := 2; N := (n-1)^n; j := 'j'; am2 := sum( (x[j]-x[i[j]])^2,j=1..n ); j := 1; k := 1; s1 := am2^(-(n-1)/2); s2 := am2^(-n/2); while k < N do am2 := the_next(i, j, n, am2, x); j := (j mod n) + 1; od; end; # exact_Eh npde[e_Eh] :=proc (x,n) local i,a,j,k,N,n1,n2,am2,s1,s2,C; for j to n do i[j] := 1; od; i[1] := 2; N := (n-1)^n; j := 'j'; am2 := sum( (x[j]-x[i[j]])^2,j=1..n ); j := 1; k := 1; n1 := (n-1)/2; n2 := n/2; s1 := am2^(-n1); s2 := am2^(-n2); for k from 2 by 1 to N do am2 := Next(i, n, am2, x); s1 := s1 + am2^(-n1); s2 := s2 + am2^(-n2); od; C := evalf( GAMMA((n-1)/2)/GAMMA(n/2)/sqrt(2) ); RETURN(evalf( C*s1/s2 )); end; # e_Eh npde[Next] :=proc (i,n,am2,x) local j,r; r := am2; j := 1; while i[j] = n and j < n do r := r + inc(i,j,n,x); j := j+1; od; if j <= n then r := r + inc(i,j,n,x); fi; RETURN(r); end; # Next npde[the_next] :=proc (i,j,n,am2,x) local oldx,newx; print(`j`=j,`n`=n); oldx := x[i[j]]; i[j] := (i[j] + 1) mod n; if i[j] = j then i[j] := (j+1) mod n; fi; print('ok2', `j`=j); newx := x[i[j]]; RETURN(am2 + (newx - oldx)*(newx+oldx - 2*oldx)); end; # the_next npde[inc] :=proc (i,j,n,x) local k,old_del,new_del; k := i[j]; old_del := (x[j]-x[k])^2; if k <> (j-1) and k <> n then i[j] := k+1; fi; if k = n and j <> 1 then i[j] := 1; fi; if k = n and j = 1 then i[j] := 2; fi; if k = (j-1) and j <> n then i[j] := j+1; fi; if k = (j-1) and j = n then i[j] := 1; fi; new_del := (x[j] - x[i[j]])^2; RETURN( evalf(new_del - old_del) ); end; # inc npde[test]:=proc () n := 5; x := sort([stats[random,normald[0,1]](n)]); i := vector([seq((j mod n) + 1,j=1..n)]); j := 'j'; am2 := sum( (x[j]-x[i[j]])^2,j=1..n ); an := am2^(-n/2); L := one_iter(an,am2,x,i,n); h_hat_exact := e_Eh(x,n); end; # npde[test] npde[test2]:=proc (n) local i,j,am2,an; global x,h_hat_exact,n_iter,n_burn,n_inter; x := sort([stats[random,normald[0,1]](n)]); i := vector([seq((j mod n) + 1,j=1..n)]); j := 'j'; am2 := sum( (x[j]-x[i[j]])^2,j=1..n ); an := am2^(-n/2); h_hat_exact := e_Eh(x,n); EhA(n_iter,n_burn,n_inter,x,n); end; # npde[test2] npde[cv_likelihood]:=proc (X::vector,n::{posint,name},h::{realcons,name}) local p,s,i,j; p := 1; s := 1; for i to n do s := 0; for j to n do if j <> i then s := s + exp(-((X[i] - X[j])/h)^2/2); fi; od; p := p*s; od; RETURN( p/h^n ); end; # npde[cv_likelihood]