Directory[]

SetDirectory["/Users/roberto/cernbox/Working/Teaching/Monte Carlo/Lezioni"]

<<FixPolygons`

Get["/Users/roberto/cernbox/Working/Teaching/Metodi/Lezioni/"<>"utility.m"]

?replaceGenericPoints

?GenericPoints

replaceGenericPoints[listofpoints_]:=Module[{repl,len=Length[listofpoints]},repl=Flatten[Table[{x[i]->listofpoints[[All,1]][[i]],y[i]->listofpoints[[All,2]][[i]]},{i,1,len}]]]

GenericPoints[P_]:=Table[{x[k],y[k]},{k,1,P}]

lossFunction[n_]:=Collect[Total[Map[(model[#[[1]]]-#[[2]])^2&,GenericPoints[n]]],{a,b},Simplify]

RealQ[x_]:=Element[x,Reals]===True

Clear[MinimizeLoss]
Options[MinimizeLoss]={Parameters->{a,b}}
MinimizeLoss[points_,opt:OptionsPattern[]]:=Module[{valsol,ssols,res,$loss,eqA,eqB,par1,par2,sols},
(*Print[OptionValue[Parameters]];*)
par1=(OptionValue[Parameters])[[1]];
par2=(OptionValue[Parameters])[[2]];
$loss=ReplaceAll[lossFunction[Length[points]],replaceGenericPoints[points]];
eqA=D[$loss,par1]==0;
eqB=D[$loss,par2]==0;
sols=NSolve[{eqA,eqB},{par1,par2}];(*{par1,par2}/.*)
(*Print[sols];*)
ssols=Select[sols,(RealQ[par1/.#]&&RealQ[par2/.#])&];
(*Print[ssols];*)
valsol=Map[{Chop[ReplaceAll[$loss,#]],#}&,ssols];
(*Print[valsol];*)
res=SortBy[valsol,First];
(*Print["res=",{par1,par2}/.res[[1,2]]];*)
{par1,par2}/.res[[1,2]]
]

Contour2D[betfitS_]:=Module[{xyz,xy,z,histoL},
histoL=HistogramList[Take[betfitS,{1,-1}]];
z=Flatten[histoL[[2]],1];
xy=Tuples[{N[ListToMidPoints[histoL[[1,1]]]],N[ListToMidPoints[histoL[[1,2]]]]}];
xyz=MapThread[Append,{xy,z}];
ListContourPlot[xyz,ContourLabels->True,(*Contours->{6700,9500,9900},*)PlotRange->All]
]

Straight Line

model[x_]:=a+b*x

Generically a point does not lie on the model curve, in our case on the straight line. So we associate a “loss” to the fact that  |f(x)-y| > 0. The choice of the loss function is our choice, depends on the application.

For usual curve - fitting of points that we presume to be random variable extracted from a probability density function close to our "model" f (x) one uses the square distance in units of the uncertainty

For just two points the equation is pretty simple

GenericPoints[2]
Map[(model[#[[1]]]-#[[2]])^2&,%]//Total
loss=Collect[%,{a,b},Simplify]

We can find a “minimum loss” point by looking for extremal points (first derivative =0) and can be solved (just one solution, we should check it’s a minimum...).

D[loss,a]==0
D[loss,b]==0
Solve[{%,%%},{a,b}]

For n-points is more complex, but not more complicated than a 2-nd order equation in each of the parameters.

GenericPoints[3]
Map[(model[#[[1]]]-#[[2]])^2&,%]//Total
loss=Collect[%,{a,b},Simplify]

D[loss,a]==0
D[loss,b]==0
Solve[{%,%%},{a,b}]

points={{1,1},{2,2},{3,4},{4,4.9}}
$loss=Collect[ReplaceAll[lossFunction[Length[points]],replaceGenericPoints[points]],{a,b}]
fit=NonlinearModelFit[points,model[x],{a,b},x];
ContourPlot[$loss,{a,-3*0.4-0.45,0.45+3*0.4},{b,-3*0.2+1.37,1.37+3*0.2},Contours->{0.5,1,3,5,7,10},ContourLabels->(Text[Framed[#3],{#1,#2},Background->White]&),FrameLabel->{"a","b"},Epilog->{Text[Style["X",Bold,Red],{a,b}/.fit["BestFitParameters"]] }]
Show[ListPlot[points],Plot[fit["BestFit"],{x,0,5}]]

fit["ParameterErrors"]

fit["EstimatedVariance"]

fit["FitResiduals"]
2Mean[%^2]

fit["BestFitParameters"]

fit["BestFit"]

D[$loss,a]==0
D[$loss,b]==0
Solve[{%,%%},{a,b}]

?model

replicas[[1]]

replicas[[2]]

ListPlot[{replicas[[2]],replicas[[3]]}]

points={{1,1},{2,2},{3,4},{4,4.9}}
Nreplicas=10000;
(*replicas=Transpose[Map[Tuples[{{#[[1]]},RandomVariate[NormalDistribution[#[[2]],0.7*#[[2]]],Nreplicas]}]&,points]];*)
replicas=Transpose[Map[Tuples[{{#[[1]]},RandomVariate[UniformDistribution[{#[[2]]-0.7*#[[2]],#[[2]]+0.7*#[[2]]}],Nreplicas]}]&,points]];
betfitS=Map[MinimizeLoss[#]&,replicas];
ListPlot[betfitS,PlotStyle->Red,Frame->True,AspectRatio->1,Axes->False,FrameLabel->{"a","b"}]

betfitS

Histogram[betfitS[[All,2]],10,Frame->True,FrameLabel->{"b"}]

Histogram[betfitS[[All,1]],10,Frame->True,FrameLabel->{"a"}]

betfitS[[All,1]]

points={{1,1},{2,2},{3,4},{4,4.9}}
Nreplicas=10000;
(*replicas=Transpose[Map[Tuples[{{#[[1]]},RandomVariate[NormalDistribution[#[[2]],0.7*#[[2]]],Nreplicas]}]&,points]];*)
replicas=Transpose[Map[Tuples[{{#[[1]]},RandomVariate[UniformDistribution[{#[[2]]-0.7*#[[2]],#[[2]]+0.7*#[[2]]}],Nreplicas]}]&,points]];
betfitS=Map[MinimizeLoss[#][[1]]&,replicas];
ListPlot[betfitS,PlotStyle->Red,Frame->True,AspectRatio->1,Axes->False]

betfitS[[1]]

Histogram3D[Take[betfitS,{1,-1}]]

Contour2D[betfitS]

Histogram2D[betfitS[[All,1]],betfitS[[All,2]]]

Manipulate[
ContourPlot[$loss,{a,-2scale,2scale},{b,-2*scale,2*scale},FrameLabel->{"a","b"},Epilog->{Text[Style["X",Bold,Red],{a,b}/.fit["BestFitParameters"]] }]
,{scale,1,11}]

Parabola Line

model[x_]:=a+b*x+c x^2

GenericPoints[3]
Map[(model[#[[1]]]-#[[2]])^2&,%]//Total
loss=Collect[%,{a,b},Simplify]

D[loss,a]==0
D[loss,b]==0
D[loss,c]==0
Solve[{%,%%,%%%},{a,b,c}]

points={{1,1},{2,2},{3,4},{4,4.9}}
$loss=Collect[ReplaceAll[lossFunction[Length[points]],replaceGenericPoints[points]],{a,b}]
fit=NonlinearModelFit[points,model[x],{a,b,c},x];

fit["BestFit"]

Show[ListPlot[points],Plot[fit["BestFit"],{x,0,5}]]

D[$loss,a]==0
D[$loss,b]==0
D[$loss,c]==0
Solve[{%,%%,%%%},{a,b,c}]

c=c/.fit["BestFitParameters"]
ContourPlot[$loss,{a,-2,2},{b,-2,2},FrameLabel->{"a","b"},Epilog->{Text[Style["X",Bold,Red],{a,b}/.fit["BestFitParameters"]] }]
Unset[c];

a=a/.fit["BestFitParameters"];
ContourPlot[$loss,{c,-2,2},{b,-2,2},FrameLabel->{"c","b"},Epilog->{Text[Style["X",Bold,Red],{c,b}/.fit["BestFitParameters"]] }]
Unset[a];

Oscillating around a straight line

model[x_]:=a*x+b*Sin[c*x]

Show[ListPlot[points],Plot[fit["BestFit"],{x,0,5}],Plot[model[x]/.{a->1.25,b->0.1,c->2},{x,1,5},PlotStyle->Dashed]]

GenericPoints[3]
Map[(model[#[[1]]]-#[[2]])^2&,%]//Total
loss=Collect[%,{a,b},Simplify]

D[loss,a]==0
D[loss,b]==0
D[loss,c]==0
Solve[{%,%%,%%%},{a,b,c}]

$Aborted, in General the solution may be hard to find analytically.

?lossFunction

points={{1,1},{2,2},{3,4},{4,4.9}}
$loss=Collect[ReplaceAll[lossFunction[Length[points]],replaceGenericPoints[points]],{a,b}]
fit=NonlinearModelFit[points,model[x],{a,b,c},x];

fit["BestFit"]

Show[ListPlot[points],Plot[fit["BestFit"],{x,0,5}]]

fit["EstimatedVariance"]

$loss

c=c/.fit["BestFitParameters"]
ContourPlot[$loss,{a,-2,2+2},{b,-2,2},FrameLabel->{"a","b"},Contours->{1,2,3,10,20,100,500},Epilog->{Text[Style["X",Bold,Red],{a,b}/.fit["BestFitParameters"]] }]
Unset[c];

a=a/.fit["BestFitParameters"];
ContourPlot[$loss,{c,-2,2},{b,-2,2},ContourLabels->True,Contours->{0.5,1,2,3,4,5},FrameLabel->{"c","b"},Epilog->{Text[Style["X",Bold,Red],{c,b}/.fit["BestFitParameters"]] }]
Unset[a];

Funny shapes ...

General question of how to explore this landscape of non-linear constraints

fit["BestFitParameters"]

?RealQ

D[$loss,a]==0
D[$loss,b]==0
D[$loss,c]==0
sols=NSolve[{%,%%,%%%},{a,b,c}]
Select[sols,(RealQ[a/.#]&&RealQ[b/.#])&]
ParallelMap[{Chop[ReplaceAll[$loss,#]],#}&,%]
SortBy[%,First]

A number of extremal points emerge. Can be  sorted, but more than one have same loss value.

Comparing Models

points={{1,1},{2,2},{4,4.9}}
$loss=Collect[ReplaceAll[lossFunction[Length[points]],replaceGenericPoints[points]],{a,b}]
fit=NonlinearModelFit[points,model[x],{a,b,c},x,MaxIterations->1000];

fit["BestFit"]

Show[ListPlot[points],Plot[fit["BestFit"],{x,0,5}]]

If I want to test the choice of  model[x] this is cheating/useless, because in general I will find a good solution with three points and three parameters for any model[x].

model[x_]:=a+b*x+c x^2
points={{1,1},{2,2},{4,4.9}}
$loss=Collect[ReplaceAll[lossFunction[Length[points]],replaceGenericPoints[points]],{a,b}]
fit=NonlinearModelFit[points,model[x],{a,b,c},x,MaxIterations->1000];
Show[ListPlot[points],Plot[fit["BestFit"],{x,0,5}]]
fit["BestFit"]

The two models will differ outside the domain in which I had my list of points. With more points I can distinguish which model is better suited to describe (fit!) the data.

Plot[{fit["BestFit"],Out[393]},{x,0,10}]

However, if I insist on using this particular model[x], even with (Subscript[N, points]) = (Subscript[N, parameters]) I can determine all the parameters of the model (a,b,c in this case).

“Confidence zones” with replicas of the data

Unset[MinimizeLoss]

c=c/.fit["BestFitParameters"]
MinimizeLoss[points]
Unset[c];

A gaussian distribution of the y value of the points, all  with σ=1

replicas

betfitS

c=c/.fit["BestFitParameters"]
points={{1,1},{2,2},{3,4},{4,4.9}}
Nreplicas=3000;
(*replicas=Transpose[Map[Tuples[{{#[[1]]},RandomVariate[UniformDistribution[{#[[2]]-0.7*#[[2]],#[[2]]+0.7*#[[2]]}],Nreplicas]}]&,points]];*)
replicas=Transpose[Map[Tuples[{{#[[1]]},RandomVariate[NormalDistribution[#[[2]],1],Nreplicas]}]&,points]];
betfitS=Map[MinimizeLoss[#,Parameters->{a,b}]&,replicas];
Correlation[betfitS[[All,1]],betfitS[[All,2]]]
ListPlot[betfitS,PlotStyle->Red,Frame->True,AspectRatio->1,Axes->False,PlotRange->{{-2,2+2},{-2,2}}]
Unset[c];

Correlation[Drop[betfitS[[All,1]],1],Drop[betfitS[[All,2]],-1]]

Correlation[Drop[betfitS[[All,1]],1],Drop[betfitS[[All,2]],-1]]

Still σ=1, now with flat distribution on the uncertainties ->  a similar, but clearly distinct, shape.

c=c/.fit["BestFitParameters"]
points={{1,1},{2,2},{3,4},{4,4.9}}
Nreplicas=10000;
replicas=Transpose[Map[Tuples[{{#[[1]]},RandomVariate[UniformDistribution[{#[[2]]-0.5,#[[2]]+0.5}],Nreplicas]}]&,points]];
(*replicas=Transpose[Map[Tuples[{{#[[1]]},RandomVariate[NormalDistribution[#[[2]],1],Nreplicas]}]&,points]];*)
betfitS=Map[MinimizeLoss[#,Parameters->{a,b}]&,replicas];
Correlation[betfitS[[All,1]],betfitS[[All,2]]]
ListPlot[betfitS,PlotStyle->Red,Frame->True,AspectRatio->1,Axes->False,PlotRange->{{-2,2+2},{-2,2}}]
Unset[c];

Now with relative uncertainty kept constant: σ/μ=0.1. The zone also changes inclination!

c=c/.fit["BestFitParameters"]
points={{1,1},{2,2},{3,4},{4,4.9}}
Nreplicas=10000;
(*replicas=Transpose[Map[Tuples[{{#[[1]]},RandomVariate[UniformDistribution[{#[[2]]-0.7*#[[2]],#[[2]]+0.7*#[[2]]}],Nreplicas]}]&,points]];*)
replicas=Transpose[Map[Tuples[{{#[[1]]},RandomVariate[NormalDistribution[#[[2]],0.1*#[[2]]],Nreplicas]}]&,points]];
betfitS=Map[MinimizeLoss[#,Parameters->{a,b}]&,replicas];
Correlation[betfitS[[All,1]],betfitS[[All,2]]]
ListPlot[betfitS,PlotStyle->Red,Frame->True,AspectRatio->1,Axes->False,PlotRange->{{-2,2+2},{-2,2}}]
Unset[c];

a=a/.fit["BestFitParameters"]
points={{1,1},{2,2},{3,4},{4,4.9}}
Nreplicas=500;
(*replicas=Transpose[Map[Tuples[{{#[[1]]},RandomVariate[UniformDistribution[{#[[2]]-0.7*#[[2]],#[[2]]+0.7*#[[2]]}],Nreplicas]}]&,points]];*)
replicas=Transpose[Map[Tuples[{{#[[1]]},RandomVariate[NormalDistribution[#[[2]],0.1*#[[2]]],Nreplicas]}]&,points]];
betfitS=Map[MinimizeLoss[#,Parameters->{b,c}]&,replicas];
Correlation[betfitS[[All,1]],betfitS[[All,2]]]
ListPlot[betfitS,PlotStyle->Red,Frame->True,AspectRatio->1,Axes->False,PlotRange->{{-2,2+2},{-2,2}}]
Unset[a];

Correlation[betfitS[[All,1]],betfitS[[All,2]]]

a=a/.fit["BestFitParameters"]
points={{1,1},{2,2},{3,4},{4,4.9}}
Nreplicas=100;
(*replicas=Transpose[Map[Tuples[{{#[[1]]},RandomVariate[UniformDistribution[{#[[2]]-0.7*#[[2]],#[[2]]+0.7*#[[2]]}],Nreplicas]}]&,points]];*)
replicas=Transpose[Map[Tuples[{{#[[1]]},RandomVariate[NormalDistribution[#[[2]],0.1*#[[2]]],Nreplicas]}]&,points]];
betfitS=Map[MinimizeLoss[#,Parameters->{c,b}][[1]]&,replicas];
ListPlot[betfitS,PlotStyle->Red,Frame->True,AspectRatio->1,Axes->False(*PlotRange->{{-2,2+2},{-2,2}}*)]
Unset[a];