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*)

p(Si,Si-1) = |g'[Si| * 1/Sqrt[2 Pi Theta] * E^(-(g[Si] - rho*g[Si-1] - b)/(2 Theta)))^2