<<"BarCharts`";<<"Histograms`";<<"PieCharts`"
Off[General::"spell1"]
Off[Plot::"plln"]

m = -1./2.; b=8.;σ=0.3;

num=31;

x = Table[i,{i,-10,num}]

normDist=NormalDistribution[0.,σ];
rarray=RandomReal[normDist,num];
bins=Table[-1.20+(i-1)*0.20,{i,1,15}];
Histogram[rarray,HistogramRange->{-1.20,1.20},HistogramCategories->bins]

y = Table[m x[[i]]+b + rarray[[i]],{i,1,num}];

data=Table[{x[[i]],y[[i]]},{i,1,num}];

f1 = FindFit[data,c+d t,{c,d},t]

Plot[c + d t/.f1,{t,-12,num},Epilog->Prepend[Point/@data,PointSize[0.013]],AxesOrigin->{0.,0.}]

dev = Table[y[[i]]-(c+d x[[i]] /.f1),{i,1,num}];

Histogram[dev,HistogramRange->{-1.20,1.20},HistogramCategories->bins]

chisqmin=(∑_(i=1)^num)(((((dev[[i]])/σ)))^2)

delta =(∑_(i=1)^num)(1/(σ^2)) (∑_(i=1)^num)(((x[[i]])^2)/(σ^2)) - ((((∑_(i=1)^num)((x[[i]])/(σ^2))))^2);
a1 = (1/delta)((∑_(i=1)^num)(((x[[i]])^2)/(σ^2)) (∑_(i=1)^num)((y[[i]])/(σ^2)) - (∑_(i=1)^num)((x[[i]])/(σ^2)) (∑_(i=1)^num)((x[[i]] y[[i]])/(σ^2)));
b1 = (1/delta)((∑_(i=1)^num)(1/(σ^2)) (∑_(i=1)^num)((x[[i]] y[[i]])/(σ^2)) - (∑_(i=1)^num)((x[[i]])/(σ^2)) (∑_(i=1)^num)((y[[i]])/(σ^2)));
deltaa = (Sqrt[(1/delta)(∑_(i=1)^num)(((x[[i]])^2)/(σ^2))]);
deltab = (Sqrt[(1/delta)(∑_(i=1)^num)(1/(σ^2))]);
Print["a =  ",a1, "  ± ",deltaa]
Print["b = ",b1, " ± ",deltab]

Now, let's generate the (χ^2)'s as a function of a and b, and plot them.

chisqA[ax_]:= (∑_(i=1)^num)((((y[[i]] - (ax+b1*x[[i]])))^2)/(σ^2))
chisqB[bx_]:= (∑_(i=1)^num)((((y[[i]] - (a1+bx*x[[i]])))^2)/(σ^2))

Plot[chisqA[ax],{ax,a1-deltaa,a1+deltaa},AxesOrigin->{a1,chisqmin}]
Plot[chisqB[bx],{bx,b1-deltab,b1+deltab},AxesOrigin->{b1,chisqmin}]

amin=c/.f1;
bmin=d/.f1;

fa = FindRoot[chisqA[ax]==chisqmin+1.,{ax,amin}];
errora=Abs[(ax/.fa )- amin];
fb = FindRoot[chisqB[bx]==chisqmin+1.,{bx,bmin}];
errorb=Abs[(bx/.fb )- bmin];

Print["a = ",amin  ," +/- ", errora];Print["b = ",bmin  ," +/- ", errorb];

chisqArray={};

m = -1./2.; b=8.;σ=0.3;
num=31;
x = Table[i,{i,-10,num}];
normDist=NormalDistribution[0.,σ];
Do[{rarray=RandomReal[normDist,num];
y = Table[m x[[i]]+b + rarray[[i]],{i,1,num}];
data=Table[{x[[i]],y[[i]]},{i,1,num}];
f1 = FindFit[data,c+d t,{c,d},t];
dev = Table[y[[i]]-(c+d x[[i]] /.f1),{i,1,num}];
chisq=(∑_(i=1)^num)(((((dev[[i]])/σ)))^2);
 chisqArray=Append[chisqArray,chisq]},{ii,1,100}]

bins2=Table[0+(ii-1)*2.5,{ii,1,30}]
g1=Histogram[chisqArray,HistogramCategories->bins2]

probChi[x_,ν_]:=(((x^(((1/2)(ν-2)))) (E^(-(x/2))))/((2^(ν/2))  Gamma[(ν/2)]))

g2 = Plot[280*probChi[x,ν],{x,0,2 ν},PlotRange->All]/.ν->29

Show[g1,g2]

(∫_42.556^∞)probChi[chi2,29]\[DifferentialD]chi2//N