@@ -83,6 +83,13 @@ dim_m = length(m_vec)

 dim_a = length(a_vec)

 dim_b = length(b_vec)

+err_qua_m = 0

+err_qua_a = 0

+err_qua_b = 0 

+err_lin_b = 0

+err_lin_m = 0

+err_con_b = 0

+

 #Marginalization of posterior_quad over a

 int_qua_a=array(0, c(dim_m, dim_b))

 for (i in (1:dim_m)){

@@ -128,3 +135,16 @@ for (i in (1:length(b_vec))){

 int_lin_m=array(0, dim=c(dim_b))

 I=int(b_vec, norm(posterior_con))

 int_con_b=I

+

+

+

+# The integration via fit is computationally expensive compared to the simple numerical integration 

+

+

+strqua = paste("quadratic: ax^2 + mx +b with a=",toString(best_a_qua),"+-",toString(err_qua_a),",\n m=", toString(best_m_qua),"+-",toString(err_qua_m),", b=",toString(best_b_qua),"+-",toString(err_qua_b))

+strlin = paste("linear: mx +b with m=", toString(best_m_lin),"+-",toString(err_lin_m),",\n b=",toString(best_b_lin),"+-",toString(err_lin_b))

+strcon = paste("constant: b with b=",toString(best_b_con),"+-",toString(err_con_b))

+print(strqua)

+text(0,2.8,strqua)

+text(0,2.1,strlin)

+text(0,1.6,strcon)