# Demo entry 6311148

Assignment1

Submitted by anonymous on Oct 24, 2016 at 16:20
Language: Python. Code size: 2.2 kB.

```import numpy as np
import math
from scipy import linalg
b0=np.zeros((2000,2))
b1=np.zeros((2000,3))
b2=np.zeros((2000,4))
b3=np.zeros((2000,5))
#每次用500个标准正态分布的x随机数拟合回归方程 每种情况拟合出2000个回归方程 并计算出回归系数矩阵
#此处采用广义线性模型，即y=b0+b1*x+b2*x^2+b3*x^3...+bn*x^n 最高n取到带4次方项的模型
for i in range(2000):
x=np.random.randn(500).reshape(500,1)
f=x+2*np.sin(1.5*x)
y=x+2*np.sin(1.5*x)+math.sqrt(0.2)*np.random.randn(500).reshape(500,1)
b=np.linalg.lstsq(np.hstack((np.ones((500,1)),x)),y)[0]
bb=np.linalg.lstsq(np.hstack((np.ones((500,1)),x,x**2)),y)[0]
bbb=np.linalg.lstsq(np.hstack((np.ones((500,1)),x,x**2,x**3)),y)[0]
bbbb=np.linalg.lstsq(np.hstack((np.ones((500,1)),x,x**2,x**3,x**4)),y)[0]
b0[i,:]=np.array(b).reshape(1,2)
b1[i,:]=np.array(bb).reshape(1,3)
b2[i,:]=np.array(bbb).reshape(1,4)
b3[i,:]=np.array(bbbb).reshape(1,5)
#生成新的1000个标准正态分布的x随机数作为测试集 并计算不同回归方程的预测函数值矩阵
xt=np.random.randn(1000).reshape(1000,1)
ft=xt+2*np.sin(1.5*xt)
yt=ft+math.sqrt(0.2)*np.random.randn(1000).reshape(1000,1)
xt0=np.column_stack((np.ones((1000,1)),xt))
xt1=np.column_stack((np.ones((1000,1)),xt,xt**2))
xt2=np.column_stack((np.ones((1000,1)),xt,xt**2,xt**3))
xt3=np.column_stack((np.ones((1000,1)),xt,xt**2,xt**3,xt**4))
h0=np.dot(b0,xt0.T)
h1=np.dot(b1,xt1.T)
h2=np.dot(b2,xt2.T)
h3=np.dot(b3,xt3.T)
#noise=((yt-ft)**2).mean()
#计算回归方程为y=b0+b1*x的偏差与方差
h0_mean=np.mean(h0,axis=0).reshape(1000,1)
var=np.mean(map(np.var,h0.T))
bias=((ft-h0_mean)**2).mean()
result={"bias":bias,"variance":var}
print result
#计算回归方程为y=b0+b1*x+b2*x^2的偏差与方差
h1_mean=np.mean(h1,axis=0).reshape(1000,1)
var=np.mean(map(np.var,h1.T))
bias=((ft-h1_mean)**2).mean()
result={"bias":bias,"variance":var}
print result
#计算回归方程为y=b0+b1*x+b2*x^2+b3*x^3的偏差与方差
h2_mean=np.mean(h2,axis=0).reshape(1000,1)
var=np.mean(map(np.var,h2.T))
bias=((ft-h2_mean)**2).mean()
result={"bias":bias,"variance":var}
print result
#计算回归方程为y=b0+b1*x+b2*x^2+b3*x^3+b4*x^4的偏差与方差
h3_mean=np.mean(h3,axis=0).reshape(1000,1)
var=np.mean(map(np.var,h3.T))
bias=((ft-h3_mean)**2).mean()
result={"bias":bias,"variance":var}
print result
```

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.