I'm trying to adapt the script to best fit the model function convolved with a kernel before computing the chi2.
Being not a very expert I'm having trouble to implement this. I was wondering to include my convolution line (basically the interpolation of a ParallelTable with a Numerical integral in) after having defined the modfunc in the GetChisqExpr and use a new function modfunc2 to obtain the modpoints and the Chisq.

modfunc2 should look like

modfunc2:=Interpolation[ParallelTable[{x, -2.5*Log[10.,AAf/(AAf + AAf2)*1./sigpsf0^2.*Exp[-x^2./(2.*sigpsf0^2.)]NIntegrate[yExp[-y^2./(2.*sigpsf0^2.)]BesselI[0., xy/sigpsf0^2.],{y,0.,100.},WorkingPrecision -> 10., MaxRecursion -> 10., Method -> {Automatic,"SymbolicProcessing"->0.}] + AAf2/(AAf + AAf2)*1./sigpsf02^2.*Exp[-x^2./(2.*sigpsf02^2.)]NIntegrate[yExp[-y^2./(2.*sigpsf02^2.)]BesselI[0., xy/sigpsf02^2.], {y,0.,100.}, WorkingPrecision -> 10., MaxRecursion -> 10., Method -> {Automatic,"SymbolicProcessing"->0.}]]}, {x, 0.1, 30., 0.3}],InterpolationOrder ->1.];

with all parameters expept x and y variables being assigned in the same module as numerical values.

Dear Josh,

Could you please help me with an issue that I have with your MCMC code? I need to review it completely in order to learn how yo use it in my project. In your example, in order to apply some in formation about a binary system, you have calculated the true anomaly first. I have looked everywhere, but I don't know how you have obtained the model of the true anomaly. You have:

Table[{f,
Re[2 ArcTan[((1 - e) Tan[f/2])/
Sqrt[1 - e^2]] - (e Sin[f] Sqrt[1 - e^2])/(1 +
e Cos[f])]}, {f, -Pi + 2 Pi/1000., Pi - 2 Pi/1000, 2 Pi/1000.}]

But I don't know how to obtain it. Could you please guide me how to obtain the true anomaly like you have done? Because as far as I know, the true anomaly is within a transcendent relation and can not be obtained explicitly. And what is f here in your formula?

I really appreciate your help because I am stuck in the code and cannot go any further until I have found out about the model you have used. Thank you very much in advance.

I tried working this example in Mathematica:

Clear[plogexpr];
delta[x1_] :=
NDSolve[{z'[t] == - x1*z[t], z[-1] == 1}, z, {t, -1, 0}][[1, 1, 2]]
plogexpr := -(delta[x1][0] - 1)^2;
Table[plogexpr, {x1, 0, 1, 1}];
mcmc = MCMC[plogexpr, {{x1, 0, 1, Reals}}, 100]

When I run the Table, it works fine, but then MCMC complains that the diff eq is non-numerical; the problem is in the variable x1 being inside NDSolve, it seems it does not take a numerical value; if I remove the x1 from inside the NDSolve, eg

Clear[plogexpr];
delta := NDSolve[{z'[t] == - z[t], z[-1] == 1}, z, {t, -1, 0}][[1, 1, 2]]
plogexpr := -(delta[0] - x1)^2;
mcmc = MCMC[plogexpr, {{x1, 0, 1, Reals}}, 100]

then it works fine.
Any suggestion?