Fitting to a complicated implicit function

Question

I have this data:

data = {{290, 0.112089478}, {240, 0.072525194}, {220, 
0.063863975}, {200, 0.075709626}, {170, 0.076462613}, {150, 
0.070010656}, {130, 0.065390206}, {120, 0.064750065}, {100, 
0.134944754}, {91, 0.257195474}, {78, 0.470626063}, {60, 
0.45755743}, {40, 0.358068721}, {30, 0.307480383}};

and want to fit it using a function of this form:

a x ((4 Γ^2 + y^2)/y^6) 

My only variable is x, y is a function of x and the dependency is given implicitly by:

y == Sqrt[1 + β/y (1/(Exp[y/x] - 1) + 1/2 - y/2)]

The parameters that I want to fit are {a, Γ, β}

And this is what I have done so far:

fitfunc[a_?NumericQ, Γ_?NumericQ, β_?NumericQ, x_?NumericQ] := 
 y /. FindRoot[
  a x ((4 Γ^2 + y^2)/y^6) /. 
   {y == Sqrt[1 + β/y (1/(Exp[y/x] - 1) + 1/2 - y/2)]}, {y, 1.}
 ];
 FindFit[data, fitfunc[a, Γ, β, x], {a, Γ, β}, x]

Answer 1

The wildly different orders of magnitude pointed out by @robsondenke is one of the problems. Another is that x and y are highly correlated. But the biggest problem is that the model just doesn't fit the data.

First I reformulate the model so that it is a bit more stable with respect to finding the maximum likelihood estimates of the parameters. And I use NonlinearModelFit as one can extract many more goodies than FindFit with few differences in input.

data = {{290, 0.112089478}, {240, 0.072525194}, {220, 0.063863975}, 
  {200, 0.075709626}, {170, 0.076462613}, {150, 0.070010656}, 
  {130, 0.065390206}, {120, 0.064750065}, {100, 0.134944754}, 
  {91, 0.257195474}, {78, 0.470626063}, {60, 0.45755743}, 
  {40, 0.358068721}, {30, 0.307480383}};

fitfunc[b_?NumericQ, c_?NumericQ, β_?NumericQ, x_?NumericQ] := 
  x (10 b + (c/100000) y^2)/y^6 /. FindRoot[y == Sqrt[1 + β/y (1/(Exp[y/x] - 1) + 1/2 - y/2)], {y, 1}];

sol = NonlinearModelFit[data, {fitfunc[b, c, β, x], {b > 0, c > 0, β > 0}}, {{b, 8}, {c, 1.3}, {β, 15}}, x];
sol["BestFitParameters"]
(* {b -> 8.176014208713987`,c -> 7.773373668159528`,β -> 14.584722480599252`} *)

(* Original parameters *)
Solve[{a 4 Γ^2 == 10 b, a == c/100000} /. sol["BestFitParameters"], {a, Γ}][[2]]
(* {a -> 0.00007773373668159529`,Γ -> 512.7858713433279`} *)

Show[ListPlot[data],
 Plot[fitfunc[b, c, β, x] /. sol["BestFitParameters"], {x, 1, 300}]]

It's a pretty bad fit. Is there some theoretical reason why this particular function?