I want to extract the distribution parameters for a set of data. I know what the distribution is. However I find that for different sets or even subsets of this data the values of the extracted distribution using `FindDistributionParameters[Data, Distribution]` fluctuate quite a lot. 

To test this I made a small simulation:

    TestData = RandomVariate[WeibullDistribution[84.2, 366.9, -494.9], 10000];
    
    TestDistParams = {\[Alpha]W, \[Beta]W, \[Gamma]W} /. FindDistributionParameters[TestData, WeibullDistribution[\[Alpha]W, \[Beta]W, \[Gamma]W], {{\[Alpha]W, 84.2}, {\[Beta]W, 366.9}, {\[Gamma]W, -494.9}}]
    
    Show[
    		{
    			Histogram[TestData,{"Raw", Round[Sqrt[Length[TestData]]]}],
    			Plot[PDF[WeibullDistribution[TestDistParams[[1]], TestDistParams[[2]], TestDistParams[[3]]]][x],{x, -160, -115}]
    		}
    	]

What I find is that the output of `FindDistributionParameters[..]` fluctuates quite a lot, sometimes the extracted parameters are twice as large as those input into generating the test distribution -- sometimes more so. This is even case when simulating with a large sample (10000 points) and when initialising `FindDistributionParameters` with guesses exactly as that used to generate the distribution.

What I have found to be far more robust is actually just finding the values of the centre bins with the corresponding PDF value or count, and using `NonlinearModelFit` to this.

What is the best approach of accurately extracting distribution parameters from a data set, reliably and accurately? Can one use constraints in a similar way to fitting?