Creating Data (and getting it in a useable form):

Let’s write a function to define a list of random integers of a specified range and length, and then take the sum of the numbers:

totalRandomList[long_, range_]:= 
Total[RandomInteger[range,{long}]]

Checking to see if it works:

totalRandomList[100,{0,50}]

Use table to make a list of several of these random sums:

Table[totalRandomList[100,{0,50}],{6}]

Now, we will make 10,000 of these random integer sums, each of which takes 50 random number on the range of 0 to 100. Let’s call this list rawdata:

rawdata = Table[totalRandomList[100,{0,50}],{10000}];

Let’s examine the histogram of the data. Note that it appears to be a normal distribution centered about 2500, which is the number of samples in each list, 50, times the expected value of each sample, which is 50 as well. 50*50 = 2500. Note the we called this histogram “hp.”

hp= Histogram[rawdata,{5}]

Creating another random list (and not taking the sum this time):

randomList = RandomInteger[{1,7},{16}]

The function Tally counts the number of times each element appears in the list. For a large list, this should approach a uniform, or constant, distribution.

Tally[randomList]

Sort[Tally[randomList]]

We can use this syntax of Sort[Tally[ ]] to sort the previous list, rawdata, in a way that is easier to analyze because we have the counts of all of the data points in the list:

histDataRaw = Sort[Tally[rawdata]];

histDataRaw//Short

Fitting this data: The following will be a numerical mistake (bad strategy)

Knowing this looks a lot like an Gaussian, we create a probably function with which to fit the data, something like:

model = a Exp[-b (x-c)^2]

histDataRaw[[1;;10]]

We try to use the function FindFit to fit the data. It is useful to look up the syntax of the function. You give the data in the form: {x,y}, and then the model, list the paramters you want to fit, and state the variable:

theParams =FindFit[histDataRaw,model,{a,b,c},x]

The function threw an error message “encountered a gradient that is effectively zero” this means that varying the parameters doesn’t seem to change the goodness of fit for the function at all.

Let’s plot the model to see its behavior in the vicinity of the histogram:

Plot[model/.theParams,{x,1959,3050},PlotStyle-> {Thick,Red}]

It’s constant, so apparently it did not fit correctly at all.

Fitting this data: Rescaling to get a reasonable fit...

histDataRaw//Short

Let’s make a list of just the sums encountered:

theSums = Table[First[histDataRaw[[i]]],{i,1,Length[histDataRaw]}];
Short[theSums]

Find the range and average of the sums, call these rangeX and meanX:

rangeX = Max[theSums] - Min[theSums]

meanX =Mean[theSums]//N

Let’s find the maximum value of the count of a single number, call this the height of the Gaussian:

height = Max[
Table[Last[histDataRaw[[i]]],{i,1,Length[histDataRaw]}]
]

Let’s now rescale the data, see the inline comments:

histDataScaled = Table[
{(First[histDataRaw[[i]]] - meanX)/rangeX(*take each x-value, and subtract the mean, divide by the range. The data should now range from -.5 to +.5 on the x-axis*),
Last[histDataRaw[[i]]]/height(*divide each data point by the maximum y-value. The data should no range from 0 to 1 on the y-axis*)},
{i,Length[histDataRaw]}
];

histDataScaled//Dimensions

Let’s check out a scaled plot of this data:

lp =ListPlot[histDataScaled]

Define two models, our original model and a parabolic one:

model = a Exp[-b (x-c)^2]

anotherModel = a + b (c -x)^2

Let’s try to use the FindFit function again on our new, scaled data, using the same syntax as before and the original, Gaussian model:

theParams =FindFit[histDataScaled,model,{a,b,c},x]

Note that it didn’t throw an error message this time. That’s probably a good sign.

modelFit = model/.theParams

Now, we’ll fit the parabolic function as well, and see which one works better:

anotherParams =FindFit[histDataScaled,anotherModel,{a,b,c},x]

anotherModelFit = anotherModel/.anotherParams

Note that we evaluated the two models with the fitted parameters. 

Plotting the two fitted functions below:

Show[Plot[{modelFit,anotherModelFit},{x,-.6,.6}],lp]

theScaledModel = model/.theParams

Now let’s unscale the model by plugging in the appropriate value for x and multiplying by the height:

theUnscaledModel = Simplify[height theScaledModel/.{x-> (realX - meanX)/rangeX}]

lp2=ListPlot[histDataRaw]

Let’s check the fit of the unscaled model:

Show[lp2,Plot[theUnscaledModel,{realX,Min[theSums],Max[theSums]},PlotStyle->{ Thickness[0.05],Red}]]

This is a pretty good fit, so we’re happy.

Note this is very different from the scaled model:

Plot[model/.theParams,{x,1959,3050},PlotStyle-> {Thick,Red}]

Let’s make another data set and do the same thing, this time with some roughly sinusoidal data:

someData = Table[{x + RandomReal[{-.1,.1}], (1 +RandomReal[{-.15,.1}]) Sin[ 2 Pi x + RandomReal[{-.05,.05}]]},{x,0,1,.01}];

Let’s plot it:

lp = ListPlot[someData]

A high order polynomial is usually good at fitting data, especially curvy data. If we need it to be curvier, we can probably just add more terms of higher order and get a better fit. So here is our first model, a quartic polynomial. This gives 5 parameters to fit the data.

aModel = a + b x + c x^2 + d  x^3 + e x^4

Use FindFit to fit the coefficients of the polynomial:

theFit1 = FindFit[someData,aModel,{a,b,c,d,e},x]

Look at the behavior of the model:

m1 = Plot[aModel/.theFit1,{x,0,1}]

Look at the model and the data at the same time:

Show[lp,m1]

Let’s make a new model, since this data looks like it *might* be sinusoidal. Note that this model only uses 3 paramters instead of 5. This means that it must be a simpler model in some ways.

aModel2 = a Sin[ b x + c]

Fitting this model:

theFit2 = FindFit[someData,aModel2,{a,b,c},x]

The second model plotted:

m2 = Plot[aModel2/.theFit2,{x,0,1},PlotStyle-> Red]

They both fit pretty well. Which model is actually better is probably based on if you can collect more data outside of the current range to determine whether one fit deviates from fitting properly at x’s outside of the current range or if you can come up with a plausible connection to the reasons why each model would reasonably explain the data.

One thing to also consider is that if you add enough terms to an expression, you can generally get a polynomial to fit almost any set of data. Thus, this generally isn’t a good strategy to fit data, unless you have an explanation for why the data should be governed by a polynomial.

In general, less parameters are also better if you can achieve the same fit, since this means you probably chose a better function to begin with.

Dynamic[Show[lp,m1,m2]]