For cross-checking purposes, if I want to check that a NonlinearModelFit performed in Mathematica is producing the correct standard errors, I would need to produce simulated data based on the fitted equation with known parameters. How would I go about doing that?
It would be nice if next time you put the fit equation into the question. As we discussed this yesterday, I believe your question in related to the question: Why are my fitted coefficients so well-determined?
So the idea is as follows. I use R, because I don't have Mathematica. I hope you can translate the code.
First you need to define the (population) parameters of your model. These are the values, which your fit has to find later on:
## # Define the parameters: a = 11 b = 10 c = 9 d = 8
Next you need to set up some more parameters, which defines your data structure. E.g.
Sigma = 1 # this is the "error estimate" you got from the fit nSim = 400 # number of datapoints in your sample t0 = 0 # your starting time tEnd = 12 # your end time
Now you generate the fake data:
## # Generate the fake data: set.seed(3) # not important e = rnorm(nSim, mean=0, sd=Sigma) # random error. Note that R uses sigma as the second parameter. Mathematica probably uses the variance (=Sigma^2) t = seq(t0, tEnd, length=nSim) # time y = a*cos(2*pi*t) + b*sin(2*pi*t) + c*cos(4*pi*t) + d*sin(4*pi*t) # Note that I add the random error below. ## # Check the fake data: (with and without noise) plot(t, y+e) plot(t, y)
Finally, you fit the fake data and evaluate the obtained error
## # Fit your model (without intersect): df = data.frame(t=t, y=y+e) # here I add the random error lm.out = lm(y~cos(2*pi*t) + sin(2*pi*t) + cos(4*pi*t) + sin(4*pi*t) -1, df) summary(lm.out)
In order to evaluate the error of your fitted parameters, you can loop over the above stated code and check whether or not each parameter lies within the $\pm 2 \sigma$ interval 95% of the times.