There may not be "one theta to rule them all"
Compare LogLikelihood values of the results:
In[1]:= fun[
x_] := (1 - c)*(1/(E^((-a + x - c*x)^2/(2*b^2))*(b*Sqrt[2*Pi])))
In[2]:= dist =
ProbabilityDistribution[fun[x], {x, -Infinity, Infinity}];
In[3]:= data =
RandomVariate[dist /. {a -> 0, b -> 1, c -> 1/3}, 10^3,
WorkingPrecision -> 20];
In[4]:= res =
Table[Quiet@
FindDistributionParameters[data,
dist, {{a, .2}, {b, .8}, {c, start}},
ParameterEstimator -> "MaximumLikelihood",
WorkingPrecision -> 20], {start, .1, .9, .1}]
Out[4]= {{a -> -0.019932674343740414493, b -> 1.1096360920952880069,
c -> 0.26954206522212784697}, {a -> -0.018766451751580295519,
b -> 1.0447134039438880896,
c -> 0.31227976365858204554}, {a -> -0.016674487536403968415,
b -> 0.92825542430454599255,
c -> 0.38894242441693089468}, {a -> -0.015633567435080450467,
b -> 0.87030823228157333632,
c -> 0.42708825124255829112}, {a -> -0.014913446774718998486,
b -> 0.83021968892034345328,
c -> 0.45347797919258945110}, {a -> -0.014144317149532428768,
b -> 0.78740285539643202694,
c -> 0.48166370243199674812}, {a -> -0.013811227971511007229,
b -> 0.76886004700210827827,
c -> 0.49387017410287984506}, {a -> -0.015406220411774267890,
b -> 0.85765200449162658321,
c -> 0.43541966915346278263}, {a -> -0.015791913582145423886,
b -> 0.87912323840309263789, c -> 0.42128545588647653513}}
In[5]:= LogLikelihood[dist /. #, data] & /@ res
Out[5]= {-1837.0542834440804521, -1837.0542834440804521, \
-1837.0542834440804521, -1837.0542834440804521, \
-1837.0542834440804521, -1837.0542834440804521, \
-1837.0542834440804521, -1837.0542834440804521, -1837.0542834440804521}
(*compare to the "original" distribution*)
In[6]:= LogLikelihood[dist /. {a -> 0, b -> 1, c -> 1/3}, data]
Out[6]= -1837.3811539068860566
(*the original distribution may not be the best fit to the data, even \
though it was used to generate said data*)