西瓜书 5.5 编写过程(标准BP与累计BP)
话不多说先用numpy表示出数据集
Y=['色泽','根蒂','敲声','纹理','脐部','触感','密度','含糖率','好瓜与否']
D=np.array([[2,1,2,3,3,1,0.697,0.406,1],\
[3,1,1,3,3,1,0.774,0.376,1],\
[3,1,2,3,3,1,0.634,0.264,1],\
[2,1,1,3,3,1,0.608,0.318,1],\
[1,1,2,3,3,1,0.556,0.215,1],\
[2,2,2,3,2,0,0.403,0.237,1],\
[3,2,2,2,2,0,0.481,0.149,1],\
[3,2,2,3,2,1,0.437,0.211,1],\
[3,2,1,2,2,1,0.666,0.091,0],\
[2,3,3,3,1,0,0.243,0.267,0],\
[1,3,3,1,1,1,0.245,0.057,0],\
[1,1,2,1,1,0,0.343,0.099,0],\
[2,2,2,2,3,1,0.639,0.161,0],\
[1,2,1,2,3,1,0.657,0.198,0],\
[3,2,2,3,2,0,0.360,0.370,0],\
[1,1,2,1,1,1,0.593,0.042,0],\
[2,1,1,2,2,1,0.719,0.103,0]])
lamuda=0.1
v=np.ones(shape=(8,8))
b=np.ones(shape=(8))
gama=np.ones(shape=(8))
w=np.ones(shape=(2,8))
y=np.ones(shape=(2))
sita=np.ones(shape=(2))
其中属性值的数字化图方便一律使用了值而不是向量,具体对应关系如下:
色泽:浅白 1,青绿 2,乌黑 3
根蒂:蜷缩 1,稍蜷 2,硬挺 3
敲声:沉闷 1,浊响 2,清脆 3
纹理:模糊 1,稍糊 2,清晰 3
脐部:平坦 1,稍凹 2,凹陷 3
触感:硬滑 1,软粘 0
敢使用值的原因也是这些属性基本都是存在一定的强弱关系的(清晰度、蜷曲度、凹陷度等),所以放心大胆123
先编写一个初始化函数对每个阈值和连接权进行初始化,各数组命名参考书上的图5.7。使用的是8输入8隐层2输出的网络结构,那么对于好瓜应该输出(1,0),坏瓜输出(0,1)
def initial(): for i in range(8): for j in range(8): v[i,j]=random.random() gama[i]=random.random() for i in range(2): for j in range(8): w[i,j]=random.random() sita[i]=random.random() return
然后做主函数,主体和书上P104的图5.8一样,每个函数的内容后面再说,跳出条件依然先设置是循环一定的次数。
num=100 while(num): num-=1 for k in range(16):#遍历每个输入 Y(k) G(k) ee() ref(k) print(v) print(b) print(gama) print(w) print(y) print(sita)
Y(k)的作用是根据当前遍历到的第k个样本的属性值计算输出y,分两步走,先修改每个隐层神经元的输出,再修改每个输出神经元的输出,代码如下:
def Y(a): for i in range(8):#修改每个隐层输出 b[i]=sigmoidb(a,i) for i in range(2): y[i]=sigmoidy(a,i) return
sigmoidb与sigmoidy的作用就是求当前神经元的输出,累加输入与权重的乘积后减去阈值,然后返回sigmoid函数结果,注意sigmoid函数中的-次方。
def sigmoidb(a,c): s=0 for i in range(8):#遍历第a个样本每个属性值 s=s+v[i,c]*D[a,i] s=gama[c]-s return 1/(1+math.exp(s)) def sigmoidy(a,c): s=0 for i in range(8): s=s+w[c,i]*b[i] s=sita[c]-s return 1/(1+math.exp(s))
回到主函数,下一步是更新两组推导出来的辅助计算的参数g和e,为了防止与之后要加的误差计算命名冲突就把e的更新函数改成了ee
更新方式很简单,按照书上现成的公式跑一跑就行了,具体推导过程可以自己试着写一下,难度还好。这里有个问题就是y的估计值,为了方便计算在原本的数据集中加了一列
在其中好瓜添加值0,坏瓜添加值1,与原本标记相反。代码中就表现为G中的样子。
def G(a):
for i in range(2):
g[i]=y[i]*(1-y[i])*(D[a,8+i]-y[i])
return
def ee():
for i in range(8):
s=b[i]*(1-b[i])
for j in range(2):
s=s*w[j,i]*g[j]
e[i]=s
return
然后就是根据公式更新连接权、阈值:
def ref(k):
for i in range(8):
for j in range(2):
w[j,i]=w[j,i]+lamuda*g[j]*b[i]
for i in range(2):
sita[i]=sita[i]-lamuda*g[i]
for i in range(8):
gama[i]=gama[i]-lamuda*e[i]
for j in range(8):
v[i,j]=v[i,j]+lamuda*e[j]*D[k,i]
return
主体大致就这样,循环完之后输出看一下参数就行。出现溢出的话多半是sigmoid里次方符号不对,这个计算应该没有特别大或者特别小的情况。
然后为了看眼误差,写个误差验算函数:
def E(k): s=0 for i in range(2): s=s+(y[i]-D[k,8+i])*(y[i]-D[k,8+i]) s/=2 return s
加到主函数里面,写个输出
print('第',100-num,'次',k,'样本误差:',E(k))
输出结果就不放了,可以看到误差是震荡的,而且似乎始终保持在0.16以上
这时候就要用到累积BP了。为了实现累计BP,我们要稍微修改一下前面的更新过程,如何修改呢?
先来看一下5.16,我们要用这个式子来代替5.4的式子套入到5.6-5.15的推到过程中去,来得到新的参数更新估计式。
5.16其实就是对所有样例的Ek做一个求平均,那么代入到5.10可以得到,新的参数更新估计式只要带上一个求平均的过程即可。
再看5.15,可以看到需要平均的项恰好也是gj,那就照原样计算就行。
首先更改一下主函数的结构,把eh和参数更新放到i的循环中,添加一个误差判别跳出,能得到一个较好的结果,记得修改函数的参数:
initial() num=1000 avr() Ee=1 while(num): num-=1 set() for k in range(16):#遍历每个输入 Y(k) G(k) ee() ref() aa=E() if(aa>Ee): print(aa) break Ee=aa
然后对G(k)做个手脚,使其最终得到所有样本的平均值:
def G(a): for i in range(2): g[i]=g[i]+y[i]*(1-y[i])*(D[a,8+i]-y[i]) if(a==15): g[i]/=16 return
ee不用动,但为了计算△Vih,我们需要求一下所有样本属性的平均值,存放到一个数组里:
def avr(): for i in range(8): avra[i]=0 for j in range(16): avra[i]+=D[j,i] avra[i]/=16 return
同时在程序开始时将g数组置零:
def set(): for i in range(2): g[i]=0 return
修改更新函数:
def ref(): for i in range(8): for j in range(2): w[j,i]=w[j,i]+lamuda*g[j]*b[i] for i in range(2): sita[i]=sita[i]-lamuda*g[i] for i in range(8): gama[i]=gama[i]-lamuda*e[i] for j in range(8): v[i,j]=v[i,j]+lamuda*e[j]*avra[i] return
修改求误差函数:
def E(): s=0 for j in range (16): Y(j) for i in range(2): s=s+(y[i]-D[k,8+i])*(y[i]-D[k,8+i])/2 s/=16 return s
搞定,运行后可以看到误差最终停在0.12左右(但不及时跳出会稳定在0.25左右)。