cjekel / piecewise_linear_fit_py

fit piecewise linear data for a specified number of line segments
MIT License
289 stars 59 forks source link

How to use the standard erros #9

Closed TorsteinSkar closed 5 years ago

TorsteinSkar commented 5 years ago

Hi,

First of all, thanks for a great software. I'm however struggling somewhat with the standard errors. I would like to plot both the data, a piece wise linear fit and a 95 % confidence interval for the fit, but I'm not quite sure how to actually do this. Could you for instance update the fitForSpecifiedNumberOfLineSegments.py example to also plot the 95% confidence interval?

What I have been trying is se = pwlf_Fz_dy.standard_errors()

fz_soil_xHat_95 = pwlf_Fzdy_soil.predict(brkpoints_y) + se * 1.96 fz_soil_xHat_05 = pwlf_Fzdy_soil.predict(brkpoints_y) - se * 1.96 some plotting function here

but I'm not sure if this is the correct approach.

Best regards Torstein Skår

cjekel commented 5 years ago

The standard errors are for each model parameter, and represents the standard deviation if you assume that the model is the correct form for each of the model parameters (that assumption includes that the break point locations are correct).

I still need to investigate an approach of a piecewise linear standard deviation... but in the meantime you can calculated the unbiased prediction variance and use it to generate your 95% confidence interval of your model.

I include a function to calculate the sum of squares of the residuals, which is the hard part.

Then all you need to do is calculate the unbiased prediction variance. It's usually easier though to work with standard deviations instead.

# calculate the sum of the square of the residuals
ssr = my_pwlf.fit_with_breaks(my_pwlf.fit_breaks)

# calculate the unbiased standard deviation
sigma = np.sqrt(ssr / (my_pwlf.n_data - my_pwlf.n_parameters))
# sigma can be used as a prediction variance

# plot the results
plt.figure()
plt.plot(x, y, 'o')
plt.plot(xHat, yHat, '-')
plt.fill_between(xHat, yHat - (sigma*1.96), yHat + (sigma*1.96), alpha=0.1,
                 color="r")
plt.show()

https://github.com/cjekel/piecewise_linear_fit_py/blob/master/examples/fitForSpecifiedNumberOfLineSegments_standard_deviation.py

cjekel commented 5 years ago

So I've finally implemented a prediction variance that represents the uncertainty due to the lack of data. My formula is based on one from http://www2.mae.ufl.edu/haftka/vvuq/lectures/Regression-accuracy.pptx Prediction variance is the squared version of standard error. Standard error should not be confused with standard errors.

An example on how to use the prediction variance is found in https://github.com/cjekel/piecewise_linear_fit_py/blob/master/examples/prediction_variance.py

I'll upload a new version to pypi shortly. I believe this change is what you were looking for.

TorsteinSkar commented 5 years ago

Thank you very much :)