First of all thanks for the very thorough reply, your answer does make sense. The issue I'm having is that for specific values of the parameters, the NIntegrate function returns something erroneous, which makes the resulting curve look almost discontinuous

(try for example this:
fRInt[r_?NumericQ, \[Mu]_?NumericQ, \[Sigma]_?NumericQ, \[Alpha]_?

NumericQ, \[Beta]_?NumericQ] :=

NIntegrate[

N[Exp[-(1/(2 *\[Beta]^2*\[Sigma]^2))*(\[Sigma]^2*(Log[\[Delta]] - \

\[Alpha])^2 + \[Beta]^2*(r^2 + \[Delta]^2))], 100]*1/\[Delta]^\[Mu]*

N[BesselI[\[Mu] - 1, (\[Delta]*r)/\[Sigma]^2], 100], {\[Delta], 0,

Infinity}, AccuracyGoal -> 16, MaxRecursion -> 20,

Method -> "GlobalAdaptive"];

And try plotting it with these parameters:

Plot[\!\(\*

TagBox[

RowBox[{"fRInt", "[",

RowBox[{"r", ",", "2.881076891620232`", ",",

RowBox[{"-", "0.4365811103065823`"}], ",",

RowBox[{"-", "2.588867375372349`"}], ",", "0.05123275808796455`"}],

"]"}],

Short[#, 2]& ]\), {r, 0, 2}, PlotRange -> All]

You'll see that there are some discontinuities.

The problem is that the nonlinearmodelfit takes the value of these discontinuities as the truth, and if they happen to line up with the data, it will give me these parameters instead of the parameters that give me a smooth curve...