Có vẻ như vấn đề là toán học, chứ không phải mã hóa. Tôi sẽ không giúp với con trăn, vì bạn biết rõ mình đang làm gì, nhưng tôi chắc rằng mọi người sẽ đánh giá cao việc đăng lại giải pháp của bạn
Tôi có thể nghĩ ra hai phương pháp, một phương pháp thô sơ và nhanh chóng, phương pháp kia nghiêm ngặt hơn
Phương pháp 1.
- Tạo trọng số số nguyên, nhưng đảo ngược sai số [1/lỗi], nhân với một hằng số phù hợp và làm tròn thành số nguyên gần nhất.
- Tạo tập dữ liệu mới bằng cách thêm nhiều bản sao của mỗi điểm dữ liệu, tương ứng với số nguyên trên.
- Thực hiện điều chỉnh tối thiểu bình phương trên tập dữ liệu mới này.
Rõ ràng bằng cách chọn hằng số lớn phù hợp, bạn có thể lấy trọng số khá chính xác. Ưu điểm lớn là đó là một tinh chỉnh nhỏ trên mã của bạn
Đây là sự tái tạo của khóa học Stanford Stats 191 [xem https. // trang web. standford. edu/class/stats191/], sử dụng các công cụ hệ sinh thái Python, thay vì R. Đây là bài giảng "Phép đổi và Bình phương nhỏ nhất có trọng số"
Thiết lập Sổ tay Ban đầu¶
fig, ax = plt.subplots[figsize=[8, 8]] ax.plot[data['t'], data['N_t'], 'd'] ax.set_xlabel['Time'] ax.set_ylabel['Count']6 ghi lại môi trường Python và gói hiện tại,
fig, ax = plt.subplots[figsize=[8, 8]] ax.plot[data['t'], data['N_t'], 'd'] ax.set_xlabel['Time'] ax.set_ylabel['Count']7 là trình định dạng Python ưa thích của tôi
Trong 2]
%load_ext watermark
Trong 3]
%load_ext lab_black
Trong [4]
%matplotlib inline
Tất cả quá trình nhập đều có ở đây [không phải tất cả đều được sử dụng trong Notebook này]
Trong [43]
import pandas as pd import numpy as np import seaborn as sns import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from scipy import stats from statsmodels.formula.api import ols from statsmodels.formula.api import rlm from statsmodels.formula.api import mixedlm from statsmodels.formula.api import wls import statsmodels.api as sm from patsy import dmatrices from patsy.contrasts import Treatment from patsy.contrasts import Sum from patsy.contrasts import Diff import warnings
Mô hình phi tuyến¶
Chúng tôi xem xét một tập dữ liệu trong đó mô hình tuyến tính không phù hợp
Tải và Khám phá Tập dữ liệu¶
Trong [6]
data = pd.read_csv['../data/bacteria.txt', sep='\t']
Trong [7]
data.head[]
Ra[7]
tN_t0135512211231973416645142
Trong [8]
fig, ax = plt.subplots[figsize=[8, 8]] ax.plot[data['t'], data['N_t'], 'd'] ax.set_xlabel['Time'] ax.set_ylabel['Count']
Ra[8]
Text[0, 0.5, 'Count']
Điều chỉnh mô hình tuyến tính¶
Trong [9]
fig, ax = plt.subplots[figsize=[8, 8]] ax.plot[data['t'], data['N_t'], 'd'] ax.set_xlabel['Time'] ax.set_ylabel['Count']0
fig, ax = plt.subplots[figsize=[8, 8]] ax.plot[data['t'], data['N_t'], 'd'] ax.set_xlabel['Time'] ax.set_ylabel['Count']1
Ra[9]
Kết quả hồi quy OLSDep. Biến đổi. N_tR bình phương. 0. 823Mô hình. OLSAdj. bình phương R. 0. 810Phương pháp. Số bình phương nhỏ nhất F-statistic. 60. 62Ngày. Thứ ba, ngày 07 tháng 4 năm 2020Prob [F-statistic]. 3. 01e-06Thời gian. 16. 42. 55Log-Khả năng. -76. 216Không. quan sát. 15AIC. 156. Dư 4Df. 13BICK. 157. Mẫu 8Df. 1Loại hiệp phương sai. nonrobustcoefstd errtP>. t. [0. 0250. 975]Đánh chặn259. 581022. 73011. 4200. 000210. 476308. 686t-19. 46432. 500-7. 7860. 000-24. 865-14. 063toàn bộ. 14. 123Durbin-Watson. 0. 803Prob[Xe đa năng]. 0. 001Jarque-Bera [JB]. 10. 490Xiên. 1. 648Prob[JB]. 0. 00527Kurtosis. 5. 432Cond. Không. 19. 3
Cảnh báo.
[1] Lỗi chuẩn giả định rằng ma trận hiệp phương sai của các lỗi được chỉ định chính xác.
Trực quan hóa sự phù hợp của mô hình tuyến tính của chúng tôi với dữ liệu
Trong [10]
%load_ext lab_black0
Ra[10]
%load_ext lab_black1
Chẩn đoán hồi quy¶
Các ghi chú bài giảng STATS191 có bốn biểu đồ có thể được sử dụng để đánh giá chất lượng mô hình tuyến tính. Chúng tôi xác định một hàm để hiển thị những điều này, được cung cấp đối tượng RegressionResults từ một mô hình phù hợp với
fig, ax = plt.subplots[figsize=[8, 8]] ax.plot[data['t'], data['N_t'], 'd'] ax.set_xlabel['Time'] ax.set_ylabel['Count']8
Trong [77]
%load_ext lab_black2
Trong [76]
%load_ext lab_black3
Chúng tôi thấy rằng phần dư rõ ràng không được phân phối đồng đều [Residuals vs Fitted graph]
Hơn nữa, nếu chúng ta nhìn vào đồ thị ảnh hưởng và đòn bẩy của
fig, ax = plt.subplots[figsize=[8, 8]] ax.plot[data['t'], data['N_t'], 'd'] ax.set_xlabel['Time'] ax.set_ylabel['Count']9, chúng ta sẽ thấy một ngoại lệ rõ ràng
Trong [13]
%load_ext lab_black4
Trong [14]
%load_ext lab_black5
Phù hợp với số mũ¶
Ở giai đoạn này, chúng tôi quyết định rằng một mô hình tốt hơn có thể là mô hình phân rã theo hàm mũ, trong đó nhật ký của
Text[0, 0.5, 'Count']0 có thể là tuyến tính theo thời gian
Lưu ý rằng chúng ta có thể sử dụng các hàm
Text[0, 0.5, 'Count']1 bên trong công thức
fig, ax = plt.subplots[figsize=[8, 8]] ax.plot[data['t'], data['N_t'], 'd'] ax.set_xlabel['Time'] ax.set_ylabel['Count']8
Trong [15]
%load_ext lab_black6
fig, ax = plt.subplots[figsize=[8, 8]] ax.plot[data['t'], data['N_t'], 'd'] ax.set_xlabel['Time'] ax.set_ylabel['Count']1
Ra[15]
Kết quả hồi quy OLSDep. Biến đổi. np. log[N_t]R bình phương. 0. 988Mô hình. OLSAdj. bình phương R. 0. 987Phương pháp. Số bình phương nhỏ nhất F-statistic. 1104. Ngày tháng. Thứ ba, ngày 07 tháng 4 năm 2020Prob [F-statistic]. 5. 86e-14Thời gian. 16. 42. 57Log-Khả năng. 12. 896Không. quan sát. 15AIC. -21. Dư lượng 79Df. 13BICK. -20. Mẫu 38Df. 1Loại hiệp phương sai. nonrobustcoefstd errtP>. t. [0. 0250. 975]Đánh chặn5. 97320. 06099. 9220. 0005. 8446. 102t-0. 21840. 007-33. 2220. 000-0. 233-0. 204Xe vạn năng. 0. 024Durbin-Watson. 2. 659Prob[Xe đa năng]. 0. 988Jarque-Bera [JB]. 0. 204Xiên. -0. 074Prob[JB]. 0. 903Kiểu nhọn. 2. 448Điều kiện. Không. 19. 3
Cảnh báo.
[1] Lỗi chuẩn giả định rằng ma trận hiệp phương sai của các lỗi được chỉ định chính xác.
Đánh giá mô hình hàm mũ¶
Vẽ sơ đồ các dự đoán của mô hình so với thực tế theo hai cách khác nhau, chúng ta có thể thấy chúng ta có sự phù hợp tốt hơn nhiều [đồng thời giá trị R^2 gần với 1 hơn nhiều. 0]
Trong [16]
%load_ext lab_black8
Ra[16]
Text[0, 0.5, 'Count']
Trong [17]
%matplotlib inline0
%matplotlib inline1
Ra[17]
%matplotlib inline2
Đồ thị chẩn đoán cũng trông đẹp hơn [đặc biệt là Phần dư so với Đã trang bị]
Trong [18]
%matplotlib inline3
Trong 19]
%matplotlib inline4
Trong 20]
%matplotlib inline5
Phương sai không cố định¶
Bây giờ chúng ta xem xét một tập dữ liệu, với một thuật ngữ lỗi không cố định từ quan sát này sang quan sát khác [xem bên dưới để biết mô tả về tập dữ liệu]
Mô tả biến
Y. Chi tiêu giáo dục bình quân đầu người theo nhà nước X1. Thu nhập bình quân đầu người năm 1973 theo bang X2. Tỷ trọng dân số dưới 18 tuổi X3. Tỷ trọng khu vực thành thị Vùng. Các bang nằm ở khu vực nào của đất nước
Đọc và khám phá tập dữ liệu¶
Trong [21]
%matplotlib inline6
Trong [22]
%matplotlib inline7
Hết[22]
STATEYX1X2X3Region0ME235394432550811NH231457832356412VT270401132832213MA261523330584614RI30047803038711
Thật khó chịu, chúng tôi có dấu cách trong tên cột, vì vậy hãy sửa chúng
Trong [23]
%matplotlib inline8
Hết[23]
%matplotlib inline9
Trong [24]
import pandas as pd import numpy as np import seaborn as sns import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from scipy import stats from statsmodels.formula.api import ols from statsmodels.formula.api import rlm from statsmodels.formula.api import mixedlm from statsmodels.formula.api import wls import statsmodels.api as sm from patsy import dmatrices from patsy.contrasts import Treatment from patsy.contrasts import Sum from patsy.contrasts import Diff import warnings0
Ra[24]
STATEYX1X2X3Region0ME23539443255081
Trong [25]
import pandas as pd import numpy as np import seaborn as sns import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from scipy import stats from statsmodels.formula.api import ols from statsmodels.formula.api import rlm from statsmodels.formula.api import mixedlm from statsmodels.formula.api import wls import statsmodels.api as sm from patsy import dmatrices from patsy.contrasts import Treatment from patsy.contrasts import Sum from patsy.contrasts import Diff import warnings1
Hết[25]
YX1X2X3Regioncount50. 00000050. 00000050. 00000050. 00000050. 000000mean284. 6000004675. 120000325. 740000657. 8000002. 660000std61. 340136644. 50625419. 423119145. 0163961. 061574min208. 0000003448. 000000287. 000000322. 0000001. 00000025%234. 2500004137. 250000310. 750000546. 7500002. 00000050%269. 5000004706. 000000324. 500000662. 5000003. 00000075%316. 7500005054. 250000333. 000000782. 2500003. 750000max546. 0000005889. 000000386. 000000909. 0000004. 000000
Fit Linear Model¶
Trong [26]
import pandas as pd import numpy as np import seaborn as sns import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from scipy import stats from statsmodels.formula.api import ols from statsmodels.formula.api import rlm from statsmodels.formula.api import mixedlm from statsmodels.formula.api import wls import statsmodels.api as sm from patsy import dmatrices from patsy.contrasts import Treatment from patsy.contrasts import Sum from patsy.contrasts import Diff import warnings2
Hết[26]
Kết quả hồi quy OLSDep. Biến đổi. YR bình phương. 0. 591Mô hình. OLSAdj. bình phương R. 0. Phương pháp 565. Số bình phương nhỏ nhất F-statistic. 22. 19Ngày. Thứ ba, ngày 07 tháng 4 năm 2020Prob [F-statistic]. 4. 94e-09Thời gian. 16. 43. 02Khả năng đăng nhập. -253. 89Không. quan sát. 50AIC. 515. Dư 8Df. 46BIC. 523. Mô hình 4Df. 3Loại hiệp phương sai. nonrobustcoefstd errtP>. t. [0. 0250. 975]Đánh chặn-556. 5680123. 195-4. 5180. 000-804. 547-308. 589X10. 07240. 0126. 2390. 0000. 0490. 096X21. 55210. 3154. 9320. 0000. 9192. 185X3-0. 00430. 051-0. 0830. 934-0. 1080. 099Toàn diện. 0. 719Durbin-Watson. 2. 149Prob[Xe đa năng]. 0. 698Jarque-Bera [JB]. 0. 760Xiên. 0. 264Prob[JB]. 0. 684Kurtosis. 2. 705Cond. Không. 1. 03e+05
Cảnh báo.
[1] Lỗi chuẩn giả định rằng ma trận hiệp phương sai của các lỗi được chỉ định chính xác.
[2] Số điều kiện lớn, 1. 03e+05. Điều này có thể chỉ ra rằng có
đa cộng tuyến mạnh hoặc các vấn đề về số khác.
Đánh giá mô hình¶
Biểu đồ Dư so với Đã trang bị cho thấy mức tăng rõ ràng về giá trị còn lại khi giá trị được trang bị tăng [i. e. phương sai của thuật ngữ lỗi không phải là hằng số]
Trong [27]
import pandas as pd import numpy as np import seaborn as sns import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from scipy import stats from statsmodels.formula.api import ols from statsmodels.formula.api import rlm from statsmodels.formula.api import mixedlm from statsmodels.formula.api import wls import statsmodels.api as sm from patsy import dmatrices from patsy.contrasts import Treatment from patsy.contrasts import Sum from patsy.contrasts import Diff import warnings3
Hơn nữa, chúng tôi dường như có một quan sát ngoại lệ [Bang Alaska, mã AK]
Trong [28]
import pandas as pd import numpy as np import seaborn as sns import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from scipy import stats from statsmodels.formula.api import ols from statsmodels.formula.api import rlm from statsmodels.formula.api import mixedlm from statsmodels.formula.api import wls import statsmodels.api as sm from patsy import dmatrices from patsy.contrasts import Treatment from patsy.contrasts import Sum from patsy.contrasts import Diff import warnings4
Trong [29]
import pandas as pd import numpy as np import seaborn as sns import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from scipy import stats from statsmodels.formula.api import ols from statsmodels.formula.api import rlm from statsmodels.formula.api import mixedlm from statsmodels.formula.api import wls import statsmodels.api as sm from patsy import dmatrices from patsy.contrasts import Treatment from patsy.contrasts import Sum from patsy.contrasts import Diff import warnings5
Hết[29]
import pandas as pd import numpy as np import seaborn as sns import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from scipy import stats from statsmodels.formula.api import ols from statsmodels.formula.api import rlm from statsmodels.formula.api import mixedlm from statsmodels.formula.api import wls import statsmodels.api as sm from patsy import dmatrices from patsy.contrasts import Treatment from patsy.contrasts import Sum from patsy.contrasts import Diff import warnings6
Trong [30]
import pandas as pd import numpy as np import seaborn as sns import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from scipy import stats from statsmodels.formula.api import ols from statsmodels.formula.api import rlm from statsmodels.formula.api import mixedlm from statsmodels.formula.api import wls import statsmodels.api as sm from patsy import dmatrices from patsy.contrasts import Treatment from patsy.contrasts import Sum from patsy.contrasts import Diff import warnings7
Bây giờ chúng ta trực quan hóa phần dư [được chuẩn hóa] theo vùng và có vẻ như mỗi vùng có một phương sai khác nhau của thuật ngữ lỗi
Trong [78]
import pandas as pd import numpy as np import seaborn as sns import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from scipy import stats from statsmodels.formula.api import ols from statsmodels.formula.api import rlm from statsmodels.formula.api import mixedlm from statsmodels.formula.api import wls import statsmodels.api as sm from patsy import dmatrices from patsy.contrasts import Treatment from patsy.contrasts import Sum from patsy.contrasts import Diff import warnings8
Hết[78]
import pandas as pd import numpy as np import seaborn as sns import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from scipy import stats from statsmodels.formula.api import ols from statsmodels.formula.api import rlm from statsmodels.formula.api import mixedlm from statsmodels.formula.api import wls import statsmodels.api as sm from patsy import dmatrices from patsy.contrasts import Treatment from patsy.contrasts import Sum from patsy.contrasts import Diff import warnings9
Thật khó chịu, chúng tôi thấy các mã trạng thái cũng có dấu cách, vì vậy hãy sửa các mã này
Trong [34]
data = pd.read_csv['../data/bacteria.txt', sep='\t']0
Ra[34]
data = pd.read_csv['../data/bacteria.txt', sep='\t']1
Trong [35]
data = pd.read_csv['../data/bacteria.txt', sep='\t']2
Ra[35]
data = pd.read_csv['../data/bacteria.txt', sep='\t']3
Bây giờ chúng tôi làm sạch tập dữ liệu, bằng cách loại trừ tiểu bang Alaska và tạo một bản sao của DataFrame
Sau đó chúng tôi khớp lại một mô hình tuyến tính
Trong [36]
data = pd.read_csv['../data/bacteria.txt', sep='\t']4
data = pd.read_csv['../data/bacteria.txt', sep='\t']5
Ra[36]
Kết quả hồi quy OLSDep. Biến đổi. YR bình phương. 0. 497Mô hình. OLSAdj. bình phương R. 0. 463phương pháp. Số bình phương nhỏ nhất F-statistic. 14. 80Ngày. Thứ ba, ngày 07 tháng 4 năm 2020Prob [F-statistic]. 7. 65e-07Thời gian. 17. 05. 11Khả năng đăng nhập. -242. 77Không. quan sát. 49AIC. 493. Dư 5Df. 45BIC. 501. Mô hình 1Df. 3Loại hiệp phương sai. nonrobustcoefstd errtP>. t. [0. 0250. 975]Đánh chặn-277. 5773132. 423-2. 0960. 042-544. 291-10. 864X10. 04830. 0123. 9760. 0000. 0240. 073X20. 88690. 3312. 6780. 0100. 2201. 554X30. 06680. 0491. 3540. 183-0. 0330. 166Xe buýt. 0. 883Durbin-Watson. 1. 946Prob[Xe buýt đa năng]. 0. 643Jarque-Bera [JB]. 0. 932Xiên. 0. 289Prob[JB]. 0. 628Kurtosis. 2. 652Cond. Không. 1. 23e+05
Cảnh báo.
[1] Lỗi chuẩn giả định rằng ma trận hiệp phương sai của các lỗi được chỉ định chính xác.
[2] Số điều kiện lớn, 1. 23e+05. Điều này có thể chỉ ra rằng có
đa cộng tuyến mạnh hoặc các vấn đề về số khác.
Đánh giá mô hình¶
Làm sạch dữ liệu không khắc phục được phần dư phương sai non_constant của chúng tôi
Trong [39]
data = pd.read_csv['../data/bacteria.txt', sep='\t']6
Dư địa theo Vùng vẫn chênh lệch giữa các Vùng
Trong [40]
data = pd.read_csv['../data/bacteria.txt', sep='\t']7
Ra[40]
import pandas as pd import numpy as np import seaborn as sns import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from scipy import stats from statsmodels.formula.api import ols from statsmodels.formula.api import rlm from statsmodels.formula.api import mixedlm from statsmodels.formula.api import wls import statsmodels.api as sm from patsy import dmatrices from patsy.contrasts import Treatment from patsy.contrasts import Sum from patsy.contrasts import Diff import warnings9
Bình phương nhỏ nhất có trọng số¶
Chúng tôi xác định một mô hình mới, trong đó chúng tôi đã chuẩn hóa phương sai bằng trọng số với giá trị của Y [giả sử một mô hình có số hạng sai số là N[0,eps * Y]. Chúng tôi giả sử các quan sát trong mỗi Vùng có số hạng sai số phương sai không đổi
Mảng
Text[0, 0.5, 'Count']3 sẽ giữ [1/phương sai] của thuật ngữ lỗi [trong mô hình của chúng tôi]
Trong [47]
data = pd.read_csv['../data/bacteria.txt', sep='\t']9
Nhận danh sách mã Vùng duy nhất
Trong [50]
data.head[]0
Hết[50]
data.head[]1
Đối với mỗi Khu vực, hãy tính ước tính phương sai của thuật ngữ sai số cho khu vực đó và điền trọng số cho tất cả các quan sát từ khu vực đó
Trong [52]
data.head[]2
data.head[]3
Lưu ý rằng giá trị của các trọng số ở trên phù hợp với ghi chú bài giảng STATS191
Làm vừa người mẫu với quả nặng¶
Trong [59]
data.head[]4
Hết[59]
Kết quả hồi quy WLSDep. Biến đổi. YR bình phương. 0. 757Mô hình. WLSAdj. bình phương R. 0. 741Phương pháp. Số bình phương nhỏ nhất F-statistic. 46. 75Ngày. Thứ ba, ngày 07 tháng 4 năm 2020Prob [F-statistic]. 7. 10e-14Thời gian. 18. 22. 54Log-Khả năng. -231. 58Không. quan sát. 49AIC. 471. Dư 2Df. 45BIC. 478. Mẫu 7Df. 3Loại hiệp phương sai. nonrobustcoefstd errtP>. t. [0. 0250. 975]Đánh chặn-315. 517078. 178-4. 0360. 000-472. 975-158. 059X10. 06230. 0087. 9160. 0000. 0460. 078X20. 87430. 2004. 3660. 0000. 4711. 278X30. 02940. 0340. 8570. 396-0. 0400. 098Xe vạn năng. 6. 395Durbin-Watson. 2. 023Prob[Xe đa năng]. 0. 041Jarque-Bera [JB]. 2. 405Xiên. 0. 147Prob[JB]. 0. 300Kurtosis. 1. 955Điều kiện. Không. 1. 05e+05
Cảnh báo.
[1] Lỗi chuẩn giả định rằng ma trận hiệp phương sai của các lỗi được chỉ định chính xác.
[2] Số điều kiện lớn, 1. 05e+05. Điều này có thể chỉ ra rằng có
đa cộng tuyến mạnh hoặc các vấn đề về số khác.
Các kết quả trên phù hợp với kết quả STATS191 ở hai chữ số thập phân đầu tiên [tôi nghi ngờ cảnh báo
Số điều kiện lớn, 1. 05e+05. Điều này có thể chỉ ra rằng có hiện tượng đa cộng tuyến mạnh hoặc các vấn đề về số khác
nên được chú ý
Để kiểm tra trọng số, bạn có thể lấy thông qua thuộc tính
Text[0, 0.5, 'Count']4 của RegressionResults
Trong [79]
data.head[]5
Hết[79]
data.head[]6
Có vẻ như sau khi điều chỉnh mô hình Bình phương nhỏ nhất có trọng số, Phần dư Pearson là cách tốt hơn để chuẩn hóa phần dư. Biểu đồ hình hộp bên dưới cho thấy chúng tôi đã tính đến phần lớn sự khác biệt giữa các Khu vực trong mô hình của mình, nhưng không phải tất cả