Cut_Gaussian.asy (1344B)
1 import graph; 2 real normalpdf(real x, real mu, real var){ 3 return 1/sqrt(2*pi*var) * exp(-((x-mu)^2/(4*var))); 4 } 5 real pdf(real x){ 6 return normalpdf(x,0,1); 7 } 8 real normalcdf(real to, real mean, real variance){ 9 real z = (to-mean)/sqrt(2*variance); 10 real t = 1/(1+0.3275911*abs(z)); 11 real a1 = 0.254829592; 12 real a2 = -0.284496736; 13 real a3 = 1.421413741; 14 real a4 = -1.453152027; 15 real a5 = 1.061405429; 16 real erf = 1-(((((a5*t + a4)*t) + a3)*t + a2)*t+ a1)*t*exp(-z*z); 17 real sign = 1; 18 if(z < 0){ 19 sign = -1; 20 } 21 return (1/2)*(1+sign*erf); 22 } 23 24 real cutpdf(real x, real x0, real mu, real var){ 25 if(x<x0){ 26 return 0; 27 } 28 return normalpdf(x,mu,var)*1/(1-normalcdf(x0,mu,var)); 29 } 30 31 real mean=0; 32 real variance=1; 33 34 picture p; 35 size(p,5cm,5cm,IgnoreAspect); 36 draw(p, graph(new real(real x){return normalpdf(x, mean, variance);}, mean-5*variance,mean+5*variance)); 37 xaxis(p, "$t$",BottomTop, LeftTicks); 38 yaxis(p, "$f_N(t)$",LeftRight,RightTicks(trailingzero)); 39 shipout("Cut_Gaussian", p); 40 41 real[] cuts={-2,0,2}; 42 for(real cut: cuts){ 43 picture p; 44 size(p,5cm,5cm,IgnoreAspect); 45 draw(p, graph(new real(real x){return cutpdf(x, cut, mean, variance);}, mean-5*variance,mean+5*variance)); 46 xaxis(p, "$t$",BottomTop, LeftTicks); 47 yaxis(p, "$f_{N|T \geq "+string(cut)+"}(t)$",LeftRight,RightTicks(trailingzero)); 48 shipout("Cut_Gaussian_"+string(cut), p); 49 }