(* See 15.6.5 in Numerical Recipes, 3rd edition *)

(* Normal distribution *)
SeedRandom[34567];
rand=RandomVariate[NormalDistribution[0,1],{300000}];
Histogram[rand,{0.1},"PDF",Epilog->Plot[PDF[NormalDistribution[0,1],x],{x,-4,4}][[1]],Frame->True]

(* Distribution of squared of random variables that are normally distributed = that follow a (χ^2) distribution! *)
Histogram[rand^2,{0.1},"PDF",Epilog->Plot[PDF[ChiSquareDistribution[1],x],{x,0,6},PlotRange->All][[1]],Frame->True]

(* This is the CDF of the ChiSquare-Distribution with M dof *)
t1=Assuming[Δχ2>0,Refine[CDF[ChiSquareDistribution[M],Δχ2]]]//FullSimplify

PDF[ChiSquareDistribution[M],Δχ2]

(* This is the probability inside nσs of a Normal Distro *)
t2=Integrate[(1/(σx(Sqrt[2Pi])))Exp[((-x^2)/(2σx^2))],{x,-n σx,n σx}]
t2/.n->{1,2,3}//N//TableForm

(* Test the Distro *)
((Length[Select[rand,-1<#<1&]])/(Length[rand]))//N
((Length[Select[rand,-2<#<2&]])/(Length[rand]))//N
((Length[Select[rand,-3<#<3&]])/(Length[rand]))//N

(* The two should match. Solve for ((Δχ)^2) in terms of n (desired number of sigmas) and ν (number of params) *)
Solve[t1==t2,Δχ2][[1]]

(* Values of ((Δχ)^2) for various values of n (desired number of sigmas) and M (number of params) *)
Table[N[2 InverseGammaRegularized[(M/2),0,Erf[(n/(Sqrt[2]))]],5],{n,1,5},{M,1,5}]//MatrixForm