импортировать numpy как np import matplotlib.pyplot as plt import scipy.special as sp## Размер образца. n = 50## Значения предикторов. XV = np . случайный . униформа ( низкий = - 4 , высокий = 4 , размер = n ) XV . sort ()## Матрица дизайна. X = np . единицы (( n , 2 )) X [:, 1 ] = XV## Истинные коэффициенты. бета = нп . array ([ 0 , 1. ], dtype = np . float64 )## Истинные значения ответа. EY = np . точка ( X , бета )## Наблюдаемые значения ответа. Y = EY + np . случайный . нормальный ( размер = n ) * np . sqrt ( 20 )## Получить оценки коэффициентов. u , s , vt = np . linalg . svd ( X , 0 ) v = np . транспонировать ( vt ) bhat = np . точка ( v , np . точка ( np . транспонирование ( u ), Y ) / s )## Подходящие значения. Yhat = np . точка ( X , бхат )## MSE и RMSE. MSE = (( Y - EY ) ** 2 ) . Сумма () / ( п - Х . форма [ 1 ]) Š = нп . sqrt ( MSE )## Эти множители используются при построении интервалов. XtX = np . точка ( np . транспонировать ( X ), X ) V = [ np . точка ( X [ я ,:], пр . linalg . решения ( XTX , X [ я ,:])) для I в интервале ( п )] V = нп . массив ( V )## Квантиль F, используемый при построении интервала Шеффа. QF = sp . fdtri ( Х . форма [ 1 ], п - Х . форма [ 1 ], 0.95 )## Нижняя и верхняя границы полосы Шеффа. D = s * np . SQRT ( Х . форма [ 1 ] * QF * V ) LB , UB = Yhat - D , Yhat + D## Нижняя и верхняя границы точечной полосы. D = s * np . sqrt ( 2 * V ) LBP , UBP = Yhat - D , Yhat + D## Составьте сюжет. plt . clf () plt . plot ( XV , Y , 'o' , ms = 3 , color = 'gray' ) plt . hold ( True ) a = plt . plot ( XV , EY , '-' , color = 'black' ) b = plt . сюжет ( XV , LB , '-', цвет = 'красный' ) plt . сюжет ( XV , UB , '-' , цвет = 'красный' ) c = plt . plot ( XV , LBP , '-' , color = 'blue' ) plt . plot ( XV , UBP , '-' , color = 'blue' ) d = plt . сюжет ( XV, Yhat , '-' , color = 'green' ) B = plt . Легенда ( ( , д , б , с ), ( "Истина" , "Оценка" , "95 % S imultaneous СВ" , \ "95% точечно СВ" ), 'нижний левый' ) B . draw_frame ( Ложь ) plt . ylim ([ - 20 , 15 ]) plt . gca () . set_yticks ([ - 20 , - 10 , 0 , 10 , 20 ]) plt . xlim ([ - 4 , 4 ]) plt . gca () . set_xticks ([ - 4 , - 2 , 0 , 2 , 4 ]) plt . xlabel ( "X" ) plt .ylabel ( "Y" ) plt . savefig ( "regression_confidence_band.png" ) plt . savefig ( "regression_confidence_band.svg" )