3047.sagews

f(x)=taylor(erf(x),x,0,3)

var( 'p theta' )

(p, theta)

mu(h) = h*(1-2*p);
sig = sqrt(4*p*(1-p));

g(h1,h2)=f( (mu(h1)-theta)/sig )*f( (mu(h2)-theta)/sig )

p(x)=exp(-x^2/2)/sqrt(2*pi)

q(x1,x2,y1,y2)=p(x1)*p(x2)*p(y1)*p(y2)

M0=integrate(integrate(integrate(integrate( g(x1+y1,x2+y2)*q(x1,x2,y1,y2),x1,-oo,oo),x2,-oo,oo),y1,-oo,oo),y2,-oo,oo);
M1=integrate(integrate(integrate(integrate( g(x1+y1,x1+y2)*q(x1,x2,y1,y2),x1,-oo,oo),x2,-oo,oo),y1,-oo,oo),y2,-oo,oo);
M2=integrate(integrate(integrate(integrate( g(x1+y1,x1+y1)*q(x1,x2,y1,y2),x1,-oo,oo),x2,-oo,oo),y1,-oo,oo),y2,-oo,oo);

k = (M2-2*M1+M0)

show(k.simplify_full())

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{64 \, p^{6} - 192 \, p^{5} + 240 \, p^{4} + {\left(16 \, p^{4} - 32 \, p^{3} + 24 \, p^{2} - 8 \, p + 1\right)} \theta^{2} - 160 \, p^{3} + 60 \, p^{2} - 12 \, p + 1}{4 \, {\left(\pi p^{6} - 3 \, \pi p^{5} + 3 \, \pi p^{4} - \pi p^{3}\right)}}

# Note the dependence on theta -> sparse is better
plot(k(p=0.4),theta,0,1)

var('epsilon')
show(k(p=.5-epsilon).simplify_full())

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{2916 \, {\left(65536.0 \, \epsilon^{20} - 114688.0 \, \epsilon^{18} + 86016.0 \, \epsilon^{16} - 35840.0 \, \epsilon^{14} + 8960.0 \, \epsilon^{12} - 1344.0 \, \epsilon^{10} + 112.0 \, \epsilon^{8} - 4.0 \, \epsilon^{6} - {\left(-16384.0 \, \epsilon^{18} + 28672.0 \, \epsilon^{16} - 21504.0 \, \epsilon^{14} + 8960.0 \, \epsilon^{12} - 2240.0 \, \epsilon^{10} + 336.0 \, \epsilon^{8} - 28.0 \, \epsilon^{6} + \epsilon^{4}\right)} \theta^{2}\right)}}{11.390625 \, \pi + 11943936.0 \, \pi \epsilon^{20} - 29859840.0 \, \pi \epsilon^{18} + 33592320.0 \, \pi \epsilon^{16} - 22394880.0 \, \pi \epsilon^{14} + 9797760.0 \, \pi \epsilon^{12} - 2939328.0 \, \pi \epsilon^{10} + 612360.0 \, \pi \epsilon^{8} - 87480.0 \, \pi \epsilon^{6} + 8201.25 \, \pi \epsilon^{4} - 455.625 \, \pi \epsilon^{2}}

k(p=.5-epsilon).simplify_full()

-2916*(65536.0*epsilon^20 - 114688.0*epsilon^18 + 86016.0*epsilon^16 - 35840.0*epsilon^14 + 8960.0*epsilon^12 - 1344.0*epsilon^10 + 112.0*epsilon^8 - 4.0*epsilon^6 - (-16384.0*epsilon^18 + 28672.0*epsilon^16 - 21504.0*epsilon^14 + 8960.0*epsilon^12 - 2240.0*epsilon^10 + 336.0*epsilon^8 - 28.0*epsilon^6 + epsilon^4)*theta^2)/(11.390625*pi + 11943936.0*pi*epsilon^20 - 29859840.0*pi*epsilon^18 + 33592320.0*pi*epsilon^16 - 22394880.0*pi*epsilon^14 + 9797760.0*pi*epsilon^12 - 2939328.0*pi*epsilon^10 + 612360.0*pi*epsilon^8 - 87480.0*pi*epsilon^6 + 8201.25*pi*epsilon^4 - 455.625*pi*epsilon^2)

k1=2916*epsilon^4*theta^2/11.390625*pi
show(k1)

\newcommand{\Bold}[1]{\mathbf{#1}}256.000000000000 \, \pi \epsilon^{4} \theta^{2}

plot(k1(epsilon=0.1),theta,0,1)