2020/09/21

keras手寫數字辨識_AutoEncoder

MNIST AutoEncoder

Autoencoder是一種可以將資料中的重要資訊保留下來的神經網路,有點像是資料壓縮,在做資料壓縮時,會有一個 Encoder 可以壓縮資料,另外還有一個 Decoder,可以還原資料。壓縮的過程就是用更精簡的方式保存了資料。Autoencoder 跟一般資料壓縮類似,也有 Encoder和Decoder,但 Decoder 的結果,不能確保可以完全還原。

Autoencoder會試著從測試資料,自己學習出Encoder和Decoder,並盡量讓資料在壓縮後又可以還原回去。實際上最常見的應用是,對圖片進行降噪 DeNoise。Autoencoder 是一種資料壓縮/解壓縮演算法,能夠處理 (1) 特定資料 (2) 壓縮後會遺失部分資訊,也就是無法完整還原回原本的資料 (3) 可由訓練資料中自動學習壓縮的方法。特定資料的意思跟常見的語音壓縮 mp3 不同,MPEG-2 Audio Layer III (mp3),可以處理所有的語音資料,但 AutoEncoder 訓練的結果,只能處理跟訓練資料類似的語音資料。例如處理人臉圖片的 autoencoder,無法有效處理 tree 的圖片。

不管是「Encoder」還是「Decoder」都可以調整權重,如果將Encoder+Decoder的結構建立好並搭配Input當作Output的目標答案,在訓練的過程中,Autoencoder會試著找出最好的權重來使得資訊可以盡量完整還原回去,換句話說,Autoencoder可以自行找出了Encoder和Decoder。

Encoder 的效果等同於做 Dimension Reduction,Encoder會轉換原本的資料到一個新的空間,這個空間可以比原本Features描述的空間更能精簡的描述這群數據,而中間這層Layer的數值Embedding Code就是新空間裡頭的座標,有些時候我們會用這個新空間來判斷每筆資料之間的接近程度。

理論上是無法做出一個 autoencoder,其得到的壓縮效果,能夠跟類似 jpeg, mp3 這種壓縮方法一樣好,因為我們無法取得「所有」的語音/圖片資料,進行訓練。

目前 autoencoder 有兩個實用的應用:(1) data denoising 例如圖片降噪 (2) dimensionality reduction for data visulization 對於多維度離散的資料,autoencoder 能夠學習出 data projection,功能跟 PCA (Principla Compoenent Analysis) 或 t-SNE 一樣,但效果更好。(ref: 淺談降維方法中的 PCA 與 t-SNE

Simple Autoencoder

# -*- coding: utf-8 -*-
from keras.datasets import mnist
from keras.utils import np_utils
import numpy as np
np.random.seed(10)

# 讀取 mnist 資料
(x_train, y_train), (x_test, y_test) = mnist.load_data()

# 標準化
x_train = x_train.astype('float32') / 255.
x_test = x_test.astype('float32') / 255.

# 正規化資料維度,以便 Keras 處理
x_train = x_train.reshape((len(x_train), np.prod(x_train.shape[1:])))
x_test = x_test.reshape((len(x_test), np.prod(x_test.shape[1:])))

print("x_train.shape=",x_train.shape, ", x_test.shape=",x_test.shape)

# simple autoencoder
# 建立 Autoencoder Model 並使用 x_train 資料進行訓練
from keras.layers import Input, Dense
from keras.models import Model

# this is the size of our encoded representations
encoding_dim = 32  # 32 floats -> compression of factor 24.5, assuming the input is 784 floats

# this is our input placeholder
input_img = Input(shape=(784,))
# "encoded" is the encoded representation of the input
encoded = Dense(encoding_dim, activation='relu')(input_img)
# "decoded" is the lossy reconstruction of the input
decoded = Dense(784, activation='sigmoid')(encoded)

# 建立 Model 並將 loss funciton 設為 binary cross entropy
# this model maps an input to its reconstruction
autoencoder = Model(input_img, decoded)
autoencoder.compile(optimizer='adadelta', loss='binary_crossentropy')

autoencoder.fit(x_train,
                x_train,  # Label 也設為 x_train
                epochs=25,
                batch_size=128,
                shuffle=True,
                validation_data=(x_test, x_test))

##########
# 另外製作 encoder, decoder 兩個分開的 Model
# this model maps an input to its encoded representation
encoder = Model(input_img, encoded)

# create a placeholder for an encoded (32-dimensional) input
encoded_input = Input(shape=(encoding_dim,))
# retrieve the last layer of the autoencoder model
decoder_layer = autoencoder.layers[-1]
# create the decoder model
decoder = Model(encoded_input, decoder_layer(encoded_input))

# encode and decode some digits
# note that we take them from the *test* set
encoded_imgs = encoder.predict(x_test)
decoded_imgs = decoder.predict(encoded_imgs)

# use Matplotlib
import matplotlib
matplotlib.use('agg')

import matplotlib.pyplot as plt
plt.clf()
n = 10  # how many digits we will display
plt.figure(figsize=(20, 4))
for i in range(n):
    # display original
    ax = plt.subplot(2, n, i + 1)
    plt.imshow(x_test[i].reshape(28, 28))
    plt.gray()
    ax.get_xaxis().set_visible(False)
    ax.get_yaxis().set_visible(False)

    # display reconstruction
    ax = plt.subplot(2, n, i + 1 + n)
    plt.imshow(decoded_imgs[i].reshape(28, 28))
    plt.gray()
    ax.get_xaxis().set_visible(False)
    ax.get_yaxis().set_visible(False)
plt.savefig("auto.png")

訓練過程

60000/60000 [==============================] - 4s 62us/step - loss: 0.3079 - val_loss: 0.2502
Epoch 2/25
60000/60000 [==============================] - 1s 12us/step - loss: 0.2300 - val_loss: 0.2105
Epoch 3/25
60000/60000 [==============================] - 1s 12us/step - loss: 0.2009 - val_loss: 0.1898
Epoch 4/25
60000/60000 [==============================] - 1s 13us/step - loss: 0.1839 - val_loss: 0.1756
Epoch 5/25
60000/60000 [==============================] - 1s 14us/step - loss: 0.1715 - val_loss: 0.1650
Epoch 6/25
60000/60000 [==============================] - 1s 13us/step - loss: 0.1619 - val_loss: 0.1563
Epoch 7/25
60000/60000 [==============================] - 1s 12us/step - loss: 0.1540 - val_loss: 0.1490
Epoch 8/25
60000/60000 [==============================] - 1s 12us/step - loss: 0.1473 - val_loss: 0.1429
Epoch 9/25
60000/60000 [==============================] - 1s 12us/step - loss: 0.1415 - val_loss: 0.1373
Epoch 10/25
60000/60000 [==============================] - 1s 12us/step - loss: 0.1365 - val_loss: 0.1325
Epoch 11/25
60000/60000 [==============================] - 1s 12us/step - loss: 0.1320 - val_loss: 0.1282
Epoch 12/25
60000/60000 [==============================] - 1s 12us/step - loss: 0.1279 - val_loss: 0.1243
Epoch 13/25
60000/60000 [==============================] - 1s 12us/step - loss: 0.1242 - val_loss: 0.1208
Epoch 14/25
60000/60000 [==============================] - 1s 12us/step - loss: 0.1209 - val_loss: 0.1178
Epoch 15/25
60000/60000 [==============================] - 1s 12us/step - loss: 0.1180 - val_loss: 0.1151
Epoch 16/25
60000/60000 [==============================] - 1s 12us/step - loss: 0.1155 - val_loss: 0.1127
Epoch 17/25
60000/60000 [==============================] - 1s 12us/step - loss: 0.1133 - val_loss: 0.1107
Epoch 18/25
60000/60000 [==============================] - 1s 12us/step - loss: 0.1115 - val_loss: 0.1090
Epoch 19/25
60000/60000 [==============================] - 1s 12us/step - loss: 0.1098 - val_loss: 0.1075
Epoch 20/25
60000/60000 [==============================] - 1s 12us/step - loss: 0.1085 - val_loss: 0.1062
Epoch 21/25
60000/60000 [==============================] - 1s 12us/step - loss: 0.1073 - val_loss: 0.1051
Epoch 22/25
60000/60000 [==============================] - 1s 14us/step - loss: 0.1062 - val_loss: 0.1042
Epoch 23/25
60000/60000 [==============================] - 1s 13us/step - loss: 0.1053 - val_loss: 0.1033
Epoch 24/25
60000/60000 [==============================] - 1s 15us/step - loss: 0.1046 - val_loss: 0.1026
Epoch 25/25
60000/60000 [==============================] - 1s 12us/step - loss: 0.1039 - val_loss: 0.1020

這是測試的前10 筆資料,上面是原始測試圖片,下面是經過壓縮/解壓縮後,產生的圖片,因為是簡單的 autoencoder,目前的效果還不夠好。

Sparse Autoencoder

ref: Tensorflow Day17 Sparse Autoencoder

在 encoded representation 加上 sparsity constraint

原本只有限制 hidden layer 為 32 維,在 hidden representation 加上 sparsity constraint。原本所有神經元會對所有輸入資料都有反應,但我們希望神經元只對某一些訓練資料有反應,例如神經元 A 對 5 有反應,B 只對 7 有反應。讓神經元有對每一個數字都有專業工作。

可在 loss function 加上兩項,達到這個限制

  • Sparsity Regularization
  • L2 Regularization

將 activity_regularizer 增加到 Dense Layer,並將訓練次數改為 100 次(因為增加了constraint,可以訓練更多次,而不會發生 overfitting)

# this is our input placeholder
input_img = Input(shape=(784,))
# "encoded" is the encoded representation of the input
# add a Dense layer with a L1 activity regularizer
encoded = Dense(encoding_dim, activation='relu',
                activity_regularizer=regularizers.l1(10e-5))(input_img)

# "decoded" is the lossy reconstruction of the input
decoded = Dense(784, activation='sigmoid')(encoded)

Note: 實際上這個部分測試的結果,反而變更差,目前不知道原因

Epoch 100/100
60000/60000 [==============================] - 1s 12us/step - loss: 0.2612 - val_loss: 0.2603

Deep AutoEncoder

在 encoded, decoded 從原本的一層,改為 3 層

# -*- coding: utf-8 -*-
from keras.datasets import mnist
from keras.utils import np_utils
import numpy as np
np.random.seed(10)

# 讀取 mnist 資料
(x_train, y_train), (x_test, y_test) = mnist.load_data()

# 標準化
x_train = x_train.astype('float32') / 255.
x_test = x_test.astype('float32') / 255.

# 正規化資料維度,以便 Keras 處理
x_train = x_train.reshape((len(x_train), np.prod(x_train.shape[1:])))
x_test = x_test.reshape((len(x_test), np.prod(x_test.shape[1:])))

print("x_train.shape=",x_train.shape, ", x_test.shape=",x_test.shape)

# simple autoencoder
# 建立 Autoencoder Model 並使用 x_train 資料進行訓練
from keras.layers import Input, Dense
from keras.models import Model
from keras import regularizers

# this is the size of our encoded representations
encoding_dim = 32  # 32 floats -> compression of factor 24.5, assuming the input is 784 floats

input_img = Input(shape=(784,))
encoded = Dense(128, activation='relu')(input_img)
encoded = Dense(64, activation='relu')(encoded)
encoded = Dense(32, activation='relu')(encoded)

decoded = Dense(64, activation='relu')(encoded)
decoded = Dense(128, activation='relu')(decoded)
decoded = Dense(784, activation='sigmoid')(decoded)

# 建立 Model 並將 loss funciton 設為 binary cross entropy
# this model maps an input to its reconstruction
autoencoder = Model(input_img, decoded)
autoencoder.compile(optimizer='adadelta', loss='binary_crossentropy')

autoencoder.fit(x_train,
                x_train,  # Label 也設為 x_train
                epochs=100,
                batch_size=128,
                shuffle=True,
                validation_data=(x_test, x_test))


# _________________________________________________________________
# Layer (type)                 Output Shape              Param #
# =================================================================
# input_1 (InputLayer)         (None, 784)               0
# _________________________________________________________________
# dense_1 (Dense)              (None, 128)               100480
# _________________________________________________________________
# dense_2 (Dense)              (None, 64)                8256
# _________________________________________________________________
# dense_3 (Dense)              (None, 32)                2080
# _________________________________________________________________
# dense_4 (Dense)              (None, 64)                2112
# _________________________________________________________________
# dense_5 (Dense)              (None, 128)               8320
# _________________________________________________________________
# dense_6 (Dense)              (None, 784)               101136
# =================================================================
# Total params: 222,384
# Trainable params: 222,384
# Non-trainable params: 0

decoded_imgs = autoencoder.predict(x_test)

# use Matplotlib
import matplotlib
matplotlib.use('agg')

import matplotlib.pyplot as plt
plt.clf()
n = 10  # how many digits we will display
plt.figure(figsize=(20, 4))
for i in range(n):
    # display original
    ax = plt.subplot(2, n, i + 1)
    plt.imshow(x_test[i].reshape(28, 28))
    plt.gray()
    ax.get_xaxis().set_visible(False)
    ax.get_yaxis().set_visible(False)

    # display reconstruction
    ax = plt.subplot(2, n, i + 1 + n)
    plt.imshow(decoded_imgs[i].reshape(28, 28))
    plt.gray()
    ax.get_xaxis().set_visible(False)
    ax.get_yaxis().set_visible(False)
plt.savefig("auto.png")

訓練結果 loss rate 由 0.1 降到 0.09

Epoch 100/100
60000/60000 [==============================] - 2s 33us/step - loss: 0.0927 - val_loss: 0.0925

Convolutional Autoencoder

執行前要加上環境變數

export TF_FORCE_GPU_ALLOW_GROWTH=true

執行結果

Epoch 50/50
60000/60000 [==============================] - 3s 57us/step - loss: 0.1012 - val_loss: 0.0984

Image Denoise

用加上 noise 的圖片當作 input,output 為沒有 noise 的圖片,這樣進行 model 訓練

# -*- coding: utf-8 -*-
from keras.datasets import mnist
from keras.utils import np_utils
import numpy as np
np.random.seed(10)

# 讀取 mnist 資料
(x_train, y_train), (x_test, y_test) = mnist.load_data()

# 標準化
x_train = x_train.astype('float32') / 255.
x_test = x_test.astype('float32') / 255.

# 正規化資料維度,以便 Keras 處理
x_train = np.reshape(x_train, (len(x_train), 28, 28, 1))  # adapt this if using `channels_first` image data format
x_test = np.reshape(x_test, (len(x_test), 28, 28, 1))  # adapt this if using `channels_first` image data format


# 將原圖加上 noise
noise_factor = 0.5
x_train_noisy = x_train + noise_factor * np.random.normal(loc=0.0, scale=1.0, size=x_train.shape)
x_test_noisy = x_test + noise_factor * np.random.normal(loc=0.0, scale=1.0, size=x_test.shape)

x_train_noisy = np.clip(x_train_noisy, 0., 1.)
x_test_noisy = np.clip(x_test_noisy, 0., 1.)


print("x_train.shape=",x_train.shape, ", x_test.shape=",x_test.shape)

from keras.layers import Input, Dense, Conv2D, MaxPooling2D, UpSampling2D
from keras.models import Model
from keras import backend as K

input_img = Input(shape=(28, 28, 1))  # adapt this if using `channels_first` image data format

## model

x = Conv2D(32, (3, 3), activation='relu', padding='same')(input_img)
x = MaxPooling2D((2, 2), padding='same')(x)
x = Conv2D(32, (3, 3), activation='relu', padding='same')(x)
encoded = MaxPooling2D((2, 2), padding='same')(x)

# at this point the representation is (7, 7, 32)

x = Conv2D(32, (3, 3), activation='relu', padding='same')(encoded)
x = UpSampling2D((2, 2))(x)
x = Conv2D(32, (3, 3), activation='relu', padding='same')(x)
x = UpSampling2D((2, 2))(x)
decoded = Conv2D(1, (3, 3), activation='sigmoid', padding='same')(x)

autoencoder = Model(input_img, decoded)

autoencoder.compile(optimizer='adadelta', loss='binary_crossentropy')

autoencoder.fit(x_train_noisy,
                x_train,  # Label 也設為 x_train
                epochs=100,
                batch_size=128,
                shuffle=True,
                validation_data=(x_test_noisy, x_test))


decoded_imgs = autoencoder.predict(x_test_noisy)

# use Matplotlib
import matplotlib
matplotlib.use('agg')

import matplotlib.pyplot as plt
plt.clf()
n = 10  # how many digits we will display
plt.figure(figsize=(20, 4))
for i in range(n):
    # display original
    ax = plt.subplot(2, n, i + 1)
    plt.imshow(x_test_noisy[i].reshape(28, 28))
    plt.gray()
    ax.get_xaxis().set_visible(False)
    ax.get_yaxis().set_visible(False)

    # display reconstruction
    ax = plt.subplot(2, n, i + 1 + n)
    plt.imshow(decoded_imgs[i].reshape(28, 28))
    plt.gray()
    ax.get_xaxis().set_visible(False)
    ax.get_yaxis().set_visible(False)
plt.savefig("auto.png")

結果

Epoch 100/100
60000/60000 [==============================] - 3s 56us/step - loss: 0.0941 - val_loss: 0.0941

Sequence-to-sequence Autoencoder

如果輸入的資料是 sequence,而不是 vector / 2D image,如果想要有 temporal structure 的 model 要改用 LSTM

from keras.layers import Input, LSTM, RepeatVector
from keras.models import Model

inputs = Input(shape=(timesteps, input_dim))
encoded = LSTM(latent_dim)(inputs)

decoded = RepeatVector(timesteps)(encoded)
decoded = LSTM(input_dim, return_sequences=True)(decoded)

sequence_autoencoder = Model(inputs, decoded)
encoder = Model(inputs, encoded)

Variational autoencoder (VAE)

variational autoencoder 是在 encoded representation 中增加 constraints 的 autoencoder。也就是 latent variable model,就是要學習一個原始資料統計分佈模型,接下來可用此模型,產生新的資料。

ref: variational_autoencoder.py

References

Building Autoencoders in Keras

Autoencoder 簡介與應用範例

[實戰系列] 使用 Keras 搭建一個 Denoising AE 魔法陣(模型)

機器學習技法 學習筆記 (6):神經網路(Neural Network)與深度學習(Deep Learning)

實作Tensorflow (4):Autoencoder

自編碼 Autoencoder (非監督學習)

[TensorFlow] [Keras] kernelregularizer、biasregularizer 和 activity_regularizer

2020/09/14

TensorFlow張量運算

TensorFlow 跟 Keras 最大的差異是,TensorFlow 必須要自行設計張量(矩陣)運算。

計算圖 Computational Graph

TensorFlow 設計核心是計算圖 Computational Graph,可分兩個部分:建立計算圖,及執行計算圖。

  1. 建立計算圖

    透過 TensorFlow 的模組,建立計算圖。

  2. 執行計算圖

    建立計算圖後,可產生 Session 執行計算圖。Session 的作用是在用戶端與執行裝置之間建立連結,有了連結,就可在不同裝置中,執行計算圖,後續任何跟裝置間的資料傳遞,都必須透過 Session 才能進行。

import tensorflow as tf

# 建立 tensorflow 常數, 常數值為 2, 名稱為 ts_c
ts_c = tf.constant(2,name='ts_c')
# 建立 tensorflow 變數,數值為剛剛的常數 + 5, 名稱為 ts_x
ts_x = tf.Variable(ts_c+5,name='ts_x')
print(ts_c)
## Tensor("ts_c:0", shape=(), dtype=int32)
## Tensor       就是 tensorflow 張量
## shape=()     代表這是 0 維的 tensor,也就是數值
## dtype=int32  張量資料型別為 int32
print(ts_x)
## <tf.Variable 'ts_x:0' shape=() dtype=int32_ref>


#######
# 建立 Session 執行 計算圖

# 產生 session
sess=tf.Session()
# 初始化所有 tensorflow global 變數
init = tf.global_variables_initializer()
sess.run(init)

# 透過 sess.run 執行計算圖,並列印執行結果
print('ts_c=',sess.run(ts_c))
print('ts_x=',sess.run(ts_x))

# 使用 eval,顯示 tensorflow 常數
print('ts_c=',ts_c.eval(session=sess))
print('ts_x=',ts_x.eval(session=sess))

# 不需要再使用 session 時,必須用 close 關閉 session
sess.close()

執行結果

ts_c= 2
ts_x= 7
ts_c= 2
ts_x= 7

可改用 With 語法,就不需要寫 sess.close(),會自動關閉,可解決可能沒有關閉 session 的問題,發生的原因,可能是程式忘了寫,或是中途發生錯誤。

import tensorflow as tf

# 建立 tensorflow 常數, 常數值為 2, 名稱為 ts_c
ts_c = tf.constant(2,name='ts_c')
# 建立 tensorflow 變數,數值為剛剛的常數 + 5, 名稱為 ts_x
ts_x = tf.Variable(ts_c+5,name='ts_x')
print(ts_c)
## Tensor("ts_c:0", shape=(), dtype=int32)
## Tensor       就是 tensorflow 張量
## shape=()     代表這是 0 維的 tensor,也就是數值
## dtype=int32  張量資料型別為 int32
print(ts_x)
## <tf.Variable 'ts_x:0' shape=() dtype=int32_ref>


#######
# 建立 Session 執行 計算圖

with tf.Session() as sess:
    # 初始化所有 tensorflow global 變數
    init = tf.global_variables_initializer()
    sess.run(init)

    # 透過 sess.run 執行計算圖,並列印執行結果
    print('ts_c=',sess.run(ts_c))
    print('ts_x=',sess.run(ts_x))

placeholder

剛剛建立計算圖時,常數與變數都是在建立計算圖的階段,就設定好了。但如果我們希望能在執行計算圖的階段,再設定數值,就必須使用 placeholder。

import tensorflow as tf

# 建立兩個 placeholder,然後用 multiply 相乘,結果存入 area
width = tf.placeholder("int32")
height = tf.placeholder("int32")
area=tf.multiply(width,height)

#######
# 建立 Session 執行 計算圖
# 在 sess.run 傳入 feed_dict 參數 {width: 6, height: 8}
with tf.Session() as sess:
    init = tf.global_variables_initializer()
    sess.run(init)
    print('area=',sess.run(area, feed_dict={width: 6, height: 8}))
    # area= 48

deprecated -> 改為tf.compat.v1.

import tensorflow as tf

# 建立兩個 placeholder,然後用 multiply 相乘,結果存入 area
width = tf.compat.v1.placeholder("int32")
height = tf.compat.v1.placeholder("int32")
area=tf.math.multiply(width,height)

#######
# 建立 Session 執行 計算圖
# 在 sess.run 傳入 feed_dict 參數 {width: 6, height: 8}
with tf.compat.v1.Session() as sess:
    init = tf.compat.v1.global_variables_initializer()
    sess.run(init)
    print('area=',sess.run(area, feed_dict={width: 6, height: 8}))
    # area= 48

tensorflow 數值運算方法

ref: https://www.tensorflow.org/api_docs/python/tf/math

先建立計算圖,然後用 session.run 執行

常用的數值運算

tensorflow 數值運算 說明
tf.add(x, y, name=None) 加法
tf.substract(x, y, name=None) 減法
tf.multiply(x, y, name=None) 乘法
tf.divide(x, y, name=None) 除法
tf.mod(x, y, name=None) 餘數
tf.sqrt(x, name=None) 平方
tf.abs(x, name=None) 絕對值

TensorBoard

可用視覺化的方式,查看計算圖

  1. 在建立 placeholder 與 mul 時,加上 name 參數
  2. 將 TensorBoard 的資料寫入 log file
import tensorflow as tf

# 建立兩個 placeholder,然後用 multiply 相乘,結果存入 area
width = tf.compat.v1.placeholder("int32", name="width")
height = tf.compat.v1.placeholder("int32", name="height")
area=tf.math.multiply(width,height, name="area")

#######
# 建立 Session 執行 計算圖
# 在 sess.run 傳入 feed_dict 參數 {width: 6, height: 8}
with tf.compat.v1.Session() as sess:
    init = tf.compat.v1.global_variables_initializer()
    sess.run(init)
    print('area=',sess.run(area, feed_dict={width: 6, height: 8}))
    # area= 48

    # 收集所有 TensorBoard 的資料
    tf.compat.v1.summary.merge_all()
    # 寫入 log file 到 log/area 目錄中
    train_writer = tf.compat.v1.summary.FileWriter("log/area", sess.graph)
  1. 啟動 tensorboard

    > tensorboard --logdir=~/tensorflow/log/area
    TensorBoard 1.14.0 at http://b39a314348ef:6006/ (Press CTRL+C to quit)
  2. 用 browser 瀏覽該網址,點 GRAPHS

建立 1 維與2 維張量

剛剛的例子都是0維的張量,也就是數值純量,接下來是 1 維張量 -> 向量,與2為以上的張量 -> 矩陣

dim 1 or 2 tensor

import tensorflow as tf

# 透過 tf.Variables 傳入 list 用以產生 dim 1 tensor
ts_X = tf.Variable([0.4,0.2,0.4])

# 傳入 2 維的 list 產生 dim 2 tensor
ts2_X = tf.Variable([[0.4,0.2,0.4]])

# dim 2 tensor,有三筆資料,每一筆資料有 2 個數值
W = tf.Variable([[-0.5,-0.2 ],
                 [-0.3, 0.4 ],
                 [-0.5, 0.2 ]])

with tf.compat.v1.Session() as sess:
    init = tf.compat.v1.global_variables_initializer()

    print("dim 1 tensor")
    sess.run(init)
    X=sess.run(ts_X)
    print(X)
    # [0.4 0.2 0.4]
    print("shape:", X.shape)
    # shape: (3,)

    print("")
    print("dim 2 tensor")
    X2=sess.run(ts2_X)
    print(X2)
    # [[0.4 0.2 0.4]]
    print("shape:", X2.shape)
    # shape: (1, 3)

    print("")
    print("dim 2 tensor")
    X3=sess.run(W)
    print(X3)
    print("shape:", X3.shape)
    # [[-0.5 -0.2]
    #  [-0.3  0.4]
    #  [-0.5  0.2]]
    # shape: (3, 2)

矩陣基本運算

浮點運算有近似的結果

import tensorflow as tf

# matrix multiply 矩陣乘法
X = tf.Variable([[1.,1.,1.]])

W = tf.Variable([[-0.5,-0.2 ],
                 [-0.3, 0.4 ],
                 [-0.5, 0.2 ]])

XW = tf.matmul(X,W )

# sum 加法
b1 = tf.Variable([[ 0.1,0.2]])
b2 = tf.Variable([[-1.3,0.4]])

Sum = b1+b2

# Y=X*W+b
X3 = tf.Variable([[1.,1.,1.]])

W3 = tf.Variable([[-0.5,-0.2 ],
                 [-0.3, 0.4 ],
                 [-0.5, 0.2 ]])


b3 = tf.Variable([[0.1,0.2]])

XWb =tf.matmul(X3,W3)+b3

with tf.compat.v1.Session() as sess:
    init = tf.compat.v1.global_variables_initializer()

    print("matrix multiply")
    sess.run(init)
    print(sess.run(XW ))
    # [[-1.3  0.4]]

    print("")
    print('Sum:')
    print(sess.run(Sum ))
    # [[-1.1999999  0.6      ]]

    print("")
    print('XWb:')
    print(sess.run(XWb ))
    # [[-1.1999999  0.6      ]]

以矩陣運算模擬神經網路的訊息傳導

以數學公式模擬,輸出、接收神經元的運作

\(y_1 = actication(x_1*w_{11} + x_2 * w_{21} + x_3*w_{31}+b_1)\)

\(y_2 = actication(x_1*w_{12} + x_2 * w_{22} + x_3*w_{32}+b_2)\)

合併為矩陣運算

\( \begin{bmatrix} y_1 & y_2 \end{bmatrix} = activation(\begin{bmatrix} x_1 & x_2 & x_3 \end{bmatrix} * \begin{bmatrix} w_{11} & w_{12} \\ w_{21} & w_{22} \\ w_{31} & w_{32} \end{bmatrix} + \begin{bmatrix} b_1 & b_2 \end{bmatrix} ) \)

也可表示為

\(Y = activation(X*W + b)\)

\(輸出 = 激活函數(輸入*權重 + 偏差)\)

  • 輸入 X

    有三個輸入神經元 \(x_1, x_2, x_3\),接收外部輸入

  • 輸出 y

    有兩個輸出神經元 \(y_1, y_2\)

  • 權重 W (weight)

    權重模擬神經元的軸突,連接輸入與接收(輸出)神經元,負責傳送訊息,因為要完全連接輸入與接收神經元,共需要 3(輸入) * 2(輸出) = 6 個軸突

    \(w_{11}, w_{21}, w_{31}\) 負責將 \(x_1, x_2, x_3\) 傳送訊息給 \(y_1\)

    \(w_{12}, w_{22}, w_{32}\) 負責將 \(x_1, x_2, x_3\) 傳送訊息給 \(y_2\)

  • 偏差 bias

    bias 模擬突觸瘩結構,代表接收神經元被活化的程度,偏差值越高,越容易被活化並傳遞訊息

  • 激活函數 activation function

    當接收神經元 \(y_1\) 接受刺激的總和 \(x_1*w_{11} + x_2 * w_{21} + x_3*w_{31}+b_1\) ,經過激活函數的運算,大於臨界值就會傳遞給下一個神經元

import tensorflow as tf
import numpy as np


X = tf.Variable([[0.4,0.2,0.4]])

W = tf.Variable([[-0.5,-0.2 ],
                 [-0.3, 0.4 ],
                 [-0.5, 0.2 ]])

b = tf.Variable([[0.1,0.2]])

XWb =tf.matmul(X,W)+b

# using relu actication function
# y = relu ( (X * W ) + b )
y=tf.nn.relu(tf.matmul(X,W)+b)

# using sigmoid activation function
# y = sigmoid ( (X * W ) + b )
y2=tf.nn.sigmoid(tf.matmul(X,W)+b)

with tf.compat.v1.Session() as sess:
    init = tf.compat.v1.global_variables_initializer()
    sess.run(init)
    print('X*W+b =')
    print(sess.run(XWb ))
    print('y =')
    print(sess.run(y ))
    print('y2 =')
    print(sess.run(y2 ))

    # X*W+b =
    # [[-0.35999998  0.28      ]]
    # y =
    # [[0.   0.28]]
    # y2 =
    # [[0.41095957 0.5695462 ]]

深度學習模型中,會以 Back Propagation 反向傳播演算法進行訓練, 訓練前要先建立多層感知模型,必須以亂數初始化模型的權重 weight 與 bias, tensorflow 提供 tf.random.normal 產生常態分佈的亂數矩陣

import tensorflow as tf
import numpy as np


W = tf.Variable(tf.random.normal([3, 2]))
b = tf.Variable(tf.random.normal([1, 2]))
X = tf.Variable([[0.4,0.2,0.4]])

y=tf.nn.relu(tf.matmul(X,W)+b)

with tf.compat.v1.Session() as sess:
    init = tf.compat.v1.global_variables_initializer()
    sess.run(init)
    print('b:')
    print(sess.run(b ))
    print('W:')
    print(sess.run(W ))
    print('y:')
    print(sess.run(y ))

    print('')
    # 用另一種寫法,一次取得三個 tensorflow 變數
    (b2,W2,y2)=sess.run((b,W,y))
    print('b2:')
    print(b2)
    print('W2:')
    print(W2)
    print('y2:')
    print(y2)
    # b:
    # [[0.7700923  0.02076844]]
    # W:
    # [[ 0.9547881  -0.0228505 ]
    #  [ 0.36570853  0.81177294]
    #  [ 0.0829528   0.48070174]]
    # y:
    # [[1.2583303 0.3662635]]

    # b2:
    # [[0.7700923  0.02076844]]
    # W2:
    # [[ 0.9547881  -0.0228505 ]
    #  [ 0.36570853  0.81177294]
    #  [ 0.0829528   0.48070174]]
    # y2:
    # [[1.2583303 0.3662635]]

normal distribution 的亂數

import tensorflow as tf
import numpy as np
import matplotlib.pyplot as plt

ts_norm = tf.random.normal([1000])
with tf.compat.v1.Session() as session:
    norm_data=ts_norm.eval()

print(len(norm_data))
print(norm_data[:30])

# [-0.28433087  1.4285065  -0.68437344  0.9676483  -0.80954283 -0.43311018
#   1.0973732  -1.5478781  -0.6180961  -0.9083597  -1.0577513  -0.43310425
#   0.8295066   0.80313313 -0.42189175  0.9471654  -0.00253101 -0.1117873
#   0.621246   -1.3487787  -0.79825306 -0.563185    0.50175935  0.6651971
#   1.1502678   0.2756175   0.19782086 -1.2379066   0.04300968 -1.3048639 ]

plt.hist(norm_data)
plt.savefig('normal.png')

以 placeholder 傳入 X

import tensorflow as tf
import numpy as np


W = tf.Variable(tf.random.normal([3, 2]))
b = tf.Variable(tf.random.normal([1, 2]))
X = tf.compat.v1.placeholder("float", [None,3])

y=tf.nn.relu(tf.matmul(X,W)+b)

with tf.compat.v1.Session() as sess:
    init = tf.compat.v1.global_variables_initializer()
    sess.run(init)

    X_array = np.array([[0.4,0.2,0.4]])
    (b,W,X,y)=sess.run((b,W,X,y),feed_dict={X:X_array})
    print('b:')
    print(b)
    print('W:')
    print(W)
    print('X:')
    print(X)
    print('y:')
    print(y)

    # b:
    # [[0.8461168  0.24919121]]
    # W:
    # [[ 0.10001858  0.20677406]
    #  [-0.56588995  2.555638  ]
    #  [-1.5147928  -0.43572944]]
    # X:
    # [[0.4 0.2 0.4]]
    # y:
    # [[0.16702908 0.6687367 ]]

X = tf.compat.v1.placeholder("float", [None,3]) 第1維設定為 None,是因為傳入的 X 筆數佈限數量,第 2 維是每一筆的數字個數,因為每一筆有三個數字,所以設定為 3

將 X 改為 3x3 矩陣

import tensorflow as tf
import numpy as np


W = tf.Variable(tf.random.normal([3, 2]))
b = tf.Variable(tf.random.normal([1, 2]))

X = tf.compat.v1.placeholder("float", [None,3])
y=tf.nn.relu(tf.matmul(X,W)+b)


with tf.compat.v1.Session() as sess:
    init = tf.compat.v1.global_variables_initializer()
    sess.run(init)

    X_array = np.array([[0.4,0.2 ,0.4],
                        [0.3,0.4 ,0.5],
                        [0.3,-0.4,0.5]])
    (b,W,X,y)=sess.run((b,W,X,y),feed_dict={X:X_array})
    print('b:')
    print(b)
    print('W:')
    print(W)
    print('X:')
    print(X)
    print('y:')
    print(y)

    # b:
    # [[0.6340158 0.5301216]]
    # W:
    # [[ 1.1625407   0.37071773]
    #  [-0.7906474  -0.9622891 ]
    #  [ 0.30319506 -0.04197265]]
    # X:
    # [[ 0.4  0.2  0.4]
    #  [ 0.3  0.4  0.5]
    #  [ 0.3 -0.4  0.5]]
    # y:
    # [[1.0621806  0.46916184]
    #  [0.81811655 0.23543498]
    #  [1.4506345  1.0052663 ]]

layer 函數

以相同的方式,透過 layer 函數,建立多層感知器 Multilayer perception

import tensorflow as tf
import numpy as np

# layer 函數 可用來建立 多層神經網路
# output_dim: 輸出的神經元數量
# input_dim: 輸入的神經元數量
# inputs: 輸入的 2 維陣列 placeholder
# activation: 傳入 activation function
def layer(output_dim,input_dim,inputs, activation=None):
    # 以常態分佈的亂數,建立並初始化 W weight 及 bias
    W = tf.Variable(tf.random.normal([input_dim, output_dim]))
    # 產生 (1, output_dim) 的常態分佈亂數矩陣
    b = tf.Variable(tf.random.normal([1, output_dim]))

    # 矩陣運算 XWb = (inputs * W) + b
    XWb = tf.matmul(inputs, W) + b

    # activation function
    if activation is None:
        outputs = XWb
    else:
        outputs = activation(XWb)
    return outputs

# 輸入 1x4, 第1維 因為筆數不固定,設定為 None
X = tf.placeholder("float", [None,4])
# 隱藏層 1x3
# 隱藏神經元 3 個,輸入神經元 4 個,activation function 為 relu
h = layer(output_dim=3,input_dim=4,inputs=X,
        activation=tf.nn.relu)
# 輸出 1x2
# 輸出神經元 2 個,輸入神經元 3 個,傳入 h
y = layer(output_dim=2,input_dim=3,inputs=h)
with tf.compat.v1.Session() as sess:
    init = tf.compat.v1.global_variables_initializer()
    sess.run(init)

    X_array = np.array([[0.4,0.2 ,0.4,0.5]])
    (layer_X,layer_h,layer_y)= sess.run((X,h,y),feed_dict={X:X_array})
    print('input Layer X:')
    print(layer_X)
    print('hidden Layer h:')
    print(layer_h)
    print('output Layer y:')
    print(layer_y)

    # input Layer X:
    # [[0.4 0.2 0.4 0.5]]
    # hidden Layer h:
    # [[0.9163495 0.        0.       ]]
    # output Layer y:
    # [[ 0.07022524 -2.128551  ]]

跟上面的程式功能一樣,但再加上 W, b 的 debug 資訊

import tensorflow as tf
import numpy as np

# layer 函數 可用來建立 多層神經網路
# output_dim: 輸出的神經元數量
# input_dim: 輸入的神經元數量
# inputs: 輸入的 2 維陣列 placeholder
# activation: 傳入 activation function
def layer_debug(output_dim,input_dim,inputs, activation=None):
    # 以常態分佈的亂數,建立並初始化 W weight 及 bias
    W = tf.Variable(tf.random.normal([input_dim, output_dim]))
    # 產生 (1, output_dim) 的常態分佈亂數矩陣
    b = tf.Variable(tf.random.normal([1, output_dim]))

    # 矩陣運算 XWb = (inputs * W) + b
    XWb = tf.matmul(inputs, W) + b

    # activation function
    if activation is None:
        outputs = XWb
    else:
        outputs = activation(XWb)
    return outputs, W, b

# 輸入 1x4, 第1維 因為筆數不固定,設定為 None
X = tf.placeholder("float", [None,4])
# 隱藏層 1x3
# 隱藏神經元 3 個,輸入神經元 4 個,activation function 為 relu
h,W1,b1=layer_debug(output_dim=3,input_dim=4,inputs=X,
                    activation=tf.nn.relu)
# 輸出 1x2
# 輸出神經元 2 個,輸入神經元 3 個,傳入 h
y,W2,b2=layer_debug(output_dim=2,input_dim=3,inputs=h)
with tf.compat.v1.Session() as sess:
    init = tf.compat.v1.global_variables_initializer()
    sess.run(init)

    X_array = np.array([[0.4,0.2 ,0.4,0.5]])
    (layer_X,layer_h,layer_y)= sess.run((X,h,y),feed_dict={X:X_array})
    print('input Layer X:')
    print(layer_X)
    print('W1:')
    print(W1)
    print('b1:')
    print(b1)
    print('hidden Layer h:')
    print(layer_h)
    print('W2:')
    print(W2)
    print('b2:')
    print(b2)
    print('output Layer y:')
    print(layer_y)

    # input Layer X:
    # [[0.4 0.2 0.4 0.5]]
    # W1:
    # <tf.Variable 'Variable:0' shape=(4, 3) dtype=float32_ref>
    # b1:
    # <tf.Variable 'Variable_1:0' shape=(1, 3) dtype=float32_ref>
    # hidden Layer h:
    # [[0.        0.        0.5992681]]
    # W2:
    # <tf.Variable 'Variable_2:0' shape=(3, 2) dtype=float32_ref>
    # b2:
    # <tf.Variable 'Variable_3:0' shape=(1, 2) dtype=float32_ref>
    # output Layer y:
    # [[0.68112874 0.5387946 ]]

References

TensorFlow+Keras深度學習人工智慧實務應用

2020/09/07

sippts

sippts 是用 perl 開發的一組測試 SIP Protocol 的工具。

安裝

依照 [PERL] 使用CPAN安裝模組 的說明,安裝 CPAN

yum install gcc* perl-CPAN

第一次進入 CPAN Shell 要經過一些設定,根據上面文章的說明,都是按 Enter 用預設值即可

# perl -MCPAN -e shell

根據 sippts 的說明,要安裝 CPAN modules

cpan -i IO:Socket:Timeout
cpan -i NetAddr:IP
cpan -i String:HexConvert
cpan -i Net::Address::IP::Local
cpan -i DBD::SQLite

在安裝 pcap module 前,要先安裝 libpcap

yum install libpcap libpcap-devel
cpan -f -i Net:Pcap

因為用 cpan -i Net:Pcap 安裝時,會卡在 t/04-loop.t ................ 1/195 測試,所以直接用 -f 強制安裝 pcap module

下載 sippts

wget https://github.com/Pepelux/sippts/archive/v2.0.3.tar.gz
tar zxvf v2.0.3.tar.gz
cd sippts-2.0.3

直接用 perl 執行看看,發生這樣的 error

# perl sipexten.pl
Can't locate Digest/MD5.pm

安裝 perl-Digest-MD5

yum -y install perl-Digest-MD5

使用

  • Sipscan 以 multithread 掃描 SIP services,可檢測多個 IPs, port ranges,可用在 TCP/UDP
  • Sipexten 可找到 SIP server 可運作的分機號碼。
  • Sipcracker 是 remote password cracker,可測試 SIP 密碼
  • Sipinvite 可檢測外撥電話是否需要認證,如果 SIP server 沒有做好設定,就可以外撥電話。且可以測試外撥電話後,再轉接另一個電話
  • Sipsniff 是簡單的 SIP sniffer 工具,可用 SIP method type 過濾 SIP packet
  • Sipspy 是一個簡單的 SIP server,可顯示 digest auth requests, responses
  • SipDigestLeak 可檢測是否有 Sandro Gauci 找到的 SIP digest leak vulnerability
# perl sipexten.pl -h 192.168.0.1 -e 100-200 -m REGISTER
Found match: 192.168.0.1/udp - User: 100 - Require authentication
Found match: 192.168.0.1/udp - User: 101 - Require authentication
Found match: 192.168.0.1/udp - User: 102 - Require authentication
Found match: 192.168.0.1/udp - User: 103 - Require authentication
....

sippts 預設的 useragent 為

User-Agent: pplsip

References

pplsip.SIP.Scanner

sippts v2.0.3 releases: Set of tools to audit SIP based VoIP Systems

Voip packages

2020/08/31

8D問題解決法 Eight Disciplines Problem Solving

8D問題解決法Eight Disciplines Problem Solving8D)也稱為團隊導向問題解決方法8D Report,是一個處理及解決問題的方法,常用於品質管理。

8D問題解決法的目的是識別一再出現的問題,並且要矯正並消除此問題,有助於產品及製程的提昇。若條件許可時,8D問題解決法會依照問題的統計分析來產生問題的永久對策,並且用確認根本原因的方式聚焦在問題的根源。

8D問題解決法是在汽車產業、電子組裝業界及其他產業中,利用團隊方式結構性徹底解決問題時的標準作法,常被用在回覆客戶的投訴案件報告上。

8D 是八個作業程序,品管人員依照步驟,就能夠找出問題發生的原因,並將分析問題的過程提供給客戶,同時驗證解決問題的方法,並防止問題再度發生。

8D 的起源

8D最早是美國福特公司使用。

二戰期間,美國政府率先採用一種類似8D的流程——「軍事標準1520」,又稱之為「不合格品的修正行動及部署系統」。

1987年,福特汽車公司最早以文件記錄下8D法,在其一份課程手冊中,這個方法被命名為「團隊導向的問題解決法」(Team Oriented Problem Solving)。 當時,福特的動力系統部門正被一些經年累月、反覆出現的生產問題搞得焦頭爛額,因此其管理層提請福特集團提供指導課程,幫助解決難題。

8D 的應用

  • 找出不合格的產品問題的原因及解決方法
  • 處理客戶投訴的問題的原因及解決方法
  • 分析反覆發生的問題的原因及解決方法

8D 的目標

  • 提升解決問題的效率,累積經驗
  • 提升產品品質
  • 避免或減少問題反覆發生
  • 找到問題的發生原因,提出短期,中期和長期對策並採取相應行動措施
  • 跨部門組成處理小組,改進產品開發流程的品質,防止問題再次發生

8D 的工作步驟

8D是解決問題的8條基本準則或稱8個工作步驟,但在實際應用中卻有9個步驟:

  • D0:徵兆緊急反應措施

  • D1:小組成立

  • D2:定義與描述問題

  • D3:確認、實施並確認臨時對策

  • D4:確認、識別及確認根本原因及漏失點(escape points)

  • D5:針對問題或不符合規格部份,選擇及確認永久對策

  • D6:實施永久對策

  • D7:採取預防措施

  • D8:感謝團隊成員


D0:徵兆緊急反應措施

目的:主要是為了看此類問題是否需要用 8D 來解決,如果問題太小,或是不適合用8D來解決的問題(例如價格,經費等等),這一步是針對問題發生時候的緊急反應。

關鍵:判斷問題的類型、大小、範疇等等。與D3不同,D0是問題一開始發生的反應,而D3是針對產品或服務問題本身的暫時對策。

D1:小組成立

目的:成立一個小組,小組成員具備 artifact/product 的開發知識,有足夠的時間及主導權,同時應具有所要求的能解決問題和實施對策的技術素質。小組必須有一個小組長。

關鍵:小組成員要有產品的開發知識,能夠解決問題

ex: 品管部(Quality Assurance):通常是小組召集人,負責統一回答客戶的問題。

製程部(Process):負責找出製程當中哪裡可能有問題。

測試部(Testing):尋找為何無法由測試方法檢出有問題的產品。

生產部(Production):配合工程師的要求,做實驗或收集數據,以利問題的發現並協助執行解決方案。

D2:定義與描述問題

目的:運用 5W1H (Who, What, Where, When, Why, How) 描述,來向團隊說明問題是在何時、何地、發生了什麽事、嚴重程度、目前狀態、如何緊急處理。

關鍵:搜集相關資料及數據,識別問題、確定範圍,跟客戶一起確認問題以及風險等級

D3:確認、實施並確認臨時對策

目的:保證在永久對策實施前,將問題與內外部顧客隔離。當問題發生時,不論原因找到與否,都必須要先止血,所以會先採取一些必要的暫時性措施。比如說如何在客戶端幫忙篩選(Screen, Sort)出有問題的產品,或者是更換良品給客戶,讓客戶可以繼續生產或是出貨。

在製造端應該要先採取措施防止問題產品繼續發生或出到客戶的手上,例如更換機器生產、加嚴篩選、全檢、將自動改爲手動、庫存清查等等。

暫時對策決定後,應立刻交由團隊成員帶回執行,並隨時回報成效。

關鍵:找到並選擇最佳的臨時對策,進行記錄與驗證(DOE、PPM分析、控製圖等)

D4:確認、識別及確認根本原因及漏失點(escape points)

目的:用統計工具列出可以所有潛在原因,將問題說明中提到的造成偏差的一系列事件或環境或原因相互隔離測試,並確定產生問題的根本原因。

最常使用的方法是【要因分析圖(魚骨圖)】,它提醒我們有哪些線索可以尋找,就人(員)、工(製程)、(來)料、機(器)、量(測)、及環(境)等六個面向,逐步檢討找出問題可能發生的原因,仔細比較、分析問題發生前後變動的狀況,比如說人員是否變動?作業手法是否更動?廠商來料是否變更?治具是否更換?跟環境的溫度、溼度是否相關?

經驗得知,日常作業的資料收集越齊全的工廠,找到真正問題的速度就越快,例如有日常修理報表,Cpk(統計製程),良率即時通報系統...等。

D5:針對問題或不符合規格部份,選擇及確認永久對策

目的:在生產前測試方案,並對方案進行評審以確定所選的校正對策能夠解決客戶問題,同時對其它過程不會有不良影響。

關鍵:驗證並決定最佳對策,如果有需要,就要重新評估臨時對策。將對策提交管理階層,確保能執行永久對策。

D6:實施永久對策

目的:制定一個實施永久對策的計劃,確定控制方法並納入文件,以確保消除了根本原因。在生產中應用該措施時應監督其長期效果。

關鍵:執行永久對策,廢除臨時措施。利用故障的可測量屬性,確認故障已經排除

D7:採取預防措施

目的:修改現有的管理系統、操作系統、工作慣例、設計與規程以防止這一問題與所有類似問題重覆發生。

關鍵:選擇預防措施,驗證有效性並進行監控

D8:感謝團隊成員

目的:承認小組的集體努力,對小組工作進行總結並祝賀。

關鍵:有選擇的保留重要文檔,將小組心得記錄到文件,必要的物質、精神獎勵。

References

8D問題解決法wiki

8-個解決問題的步驟

問題分析與對策解決,簡介8D report方法

8D工作方法

一步解決8D報告回復之痛

工廠8D報告

8D法是什麼?詳解8D法的九步驟!

2020/08/17

判斷點是否在多邊形內的方法

給定一個由多個點的 list 產生的多邊形,判斷另一個點座標,是否有包含在該多邊形的圖形中。

方法是從給定點座標開始,往隨便一個方向射出一條射線(例如水平往右射線),看看穿過多少條邊。如果穿過偶數次,表示點在簡單多邊形外部;如果穿過奇數次,表示點在簡單多邊形內部。

不過,要另外處理,當射線穿過頂點、射線與邊重疊的狀況,也就是給定點座標,與某一條邊共線的狀況。

這兩個連結有對方法做更詳細的說明

Point in Polygon

How to check if a given point lies inside or outside a polygon?

另外有提供一個 Java 版的實作

// A Java program to check if a given point
// lies inside a given polygon
// Refer https://www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/
// for explanation of functions onSegment(),
// orientation() and doIntersect()
class GFG
{

    // Define Infinite (Using INT_MAX
    // caused overflow problems)
    static int INF = 10000;

    static class Point
    {
        int x;
        int y;

        public Point(int x, int y)
        {
            this.x = x;
            this.y = y;
        }
    };

    // Given three colinear points p, q, r,
    // the function checks if point q lies
    // on line segment 'pr'
    static boolean onSegment(Point p, Point q, Point r)
    {
        if (q.x <= Math.max(p.x, r.x) &&
            q.x >= Math.min(p.x, r.x) &&
            q.y <= Math.max(p.y, r.y) &&
            q.y >= Math.min(p.y, r.y))
        {
            return true;
        }
        return false;
    }

    // To find orientation of ordered triplet (p, q, r).
    // The function returns following values
    // 0 --> p, q and r are colinear
    // 1 --> Clockwise
    // 2 --> Counterclockwise
    static int orientation(Point p, Point q, Point r)
    {
        int val = (q.y - p.y) * (r.x - q.x)
                - (q.x - p.x) * (r.y - q.y);

        if (val == 0)
        {
            return 0; // colinear
        }
        return (val > 0) ? 1 : 2; // clock or counterclock wise
    }

    // The function that returns true if
    // line segment 'p1q1' and 'p2q2' intersect.
    static boolean doIntersect(Point p1, Point q1,
                               Point p2, Point q2)
    {
        // Find the four orientations needed for
        // general and special cases
        int o1 = orientation(p1, q1, p2);
        int o2 = orientation(p1, q1, q2);
        int o3 = orientation(p2, q2, p1);
        int o4 = orientation(p2, q2, q1);

        // General case
        if (o1 != o2 && o3 != o4)
        {
            return true;
        }

        // Special Cases
        // p1, q1 and p2 are colinear and
        // p2 lies on segment p1q1
        if (o1 == 0 && onSegment(p1, p2, q1))
        {
            return true;
        }

        // p1, q1 and p2 are colinear and
        // q2 lies on segment p1q1
        if (o2 == 0 && onSegment(p1, q2, q1))
        {
            return true;
        }

        // p2, q2 and p1 are colinear and
        // p1 lies on segment p2q2
        if (o3 == 0 && onSegment(p2, p1, q2))
        {
            return true;
        }

        // p2, q2 and q1 are colinear and
        // q1 lies on segment p2q2
        if (o4 == 0 && onSegment(p2, q1, q2))
        {
            return true;
        }

        // Doesn't fall in any of the above cases
        return false;
    }

    // Returns true if the point p lies
    // inside the polygon[] with n vertices
    static boolean isInside(Point polygon[], int n, Point p)
    {
        // There must be at least 3 vertices in polygon[]
        if (n < 3)
        {
            return false;
        }

        // Create a point for line segment from p to infinite
        Point extreme = new Point(INF, p.y);

        // Count intersections of the above line
        // with sides of polygon
        int count = 0, i = 0;
        do
        {
            int next = (i + 1) % n;

            // Check if the line segment from 'p' to
            // 'extreme' intersects with the line
            // segment from 'polygon[i]' to 'polygon[next]'
            if (doIntersect(polygon[i], polygon[next], p, extreme))
            {
                // If the point 'p' is colinear with line
                // segment 'i-next', then check if it lies
                // on segment. If it lies, return true, otherwise false
                if (orientation(polygon[i], p, polygon[next]) == 0)
                {
                    return onSegment(polygon[i], p,
                                     polygon[next]);
                }

                count++;
            }
            i = next;
        } while (i != 0);

        // Return true if count is odd, false otherwise
        return (count % 2 == 1); // Same as (count%2 == 1)
    }

    // Driver Code
    public static void main(String[] args)
    {
        Point polygon1[] = {new Point(0, 0),
                            new Point(10, 0),
                            new Point(10, 10),
                            new Point(0, 10)};
        int n = polygon1.length;
        Point p = new Point(20, 20);
        if (isInside(polygon1, n, p))
        {
            System.out.println("Yes");
        }
        else
        {
            System.out.println("No");
        }
        p = new Point(5, 5);
        if (isInside(polygon1, n, p))
        {
            System.out.println("Yes");
        }
        else
        {
            System.out.println("No");
        }
        p = new Point(-1, 10);
        n = polygon1.length;
        if (isInside(polygon1, n, p))
        {
            System.out.println("Yes");
        }
        else
        {
            System.out.println("No");
        }


        Point polygon2[] = {new Point(0, 0),
            new Point(5, 5), new Point(5, 0)};
        p = new Point(3, 3);
        n = polygon2.length;
        if (isInside(polygon2, n, p))
        {
            System.out.println("Yes");
        }
        else
        {
            System.out.println("No");
        }
        p = new Point(5, 1);
        if (isInside(polygon2, n, p))
        {
            System.out.println("Yes");
        }
        else
        {
            System.out.println("No");
        }
        p = new Point(8, 1);
        if (isInside(polygon2, n, p))
        {
            System.out.println("Yes");
        }
        else
        {
            System.out.println("No");
        }
    }
}

以下的實作,是根據 Java 版本的內容,改成 erlang 版

-module(polygon).

%% API
-export([
  test/0,
  test2/0
]).

-record(point, {
  x :: float(),
  y :: float()
}).
-type(point() :: #point{}).

-define(INF, 10000.0).

%% 檢查 Q 是否在 線段 PR 上
%% check if point Q lise on line segment(PR)
-spec on_segment(P :: point(), Q :: point(), R :: point() ) -> boolean().
on_segment(P, Q, R) ->
  #point{x = Px, y = Py} = P,
  #point{x = Qx, y = Qy} = Q,
  #point{x = Rx, y = Ry} = R,

%%  io:format("P=~p, Q=~p, R=~p~n", [P, Q, R]),
%%
%%  io:format("Qx=~p, max(Px, Rx)=~p~n", [Qx, max(Px, Rx)]),
%%  io:format("Qx=~p, min(Px, Rx)=~p~n", [Qx, min(Px, Rx)]),
%%  io:format("Qy=~p, max(Py, Ry)=~p~n", [Qy, max(Py, Ry)]),
%%  io:format("Qy=~p, min(Py, Ry)=~p~n", [Qy, min(Py, Ry)]),
  case (Qx =< max(Px, Rx)) and (Qx >= min(Px, Rx)) and( Qy =< max(Py, Ry)) and (Qy >= min(Py, Ry)) of
    true ->
      true;
    _ ->
      false
  end.

%% 查詢 P, Q, R 的順序
%% 0: 三點共線, 1: clockwise 順時鐘, 2: counterclockwise 逆時鐘
-spec orientation(P :: point(), Q :: point(), R :: point() ) -> integer().
orientation(P, Q, R) ->
  #point{x = Px, y = Py} = P,
  #point{x = Qx, y = Qy} = Q,
  #point{x = Rx, y = Ry} = R,

  Val = (Qy - Py) * (Rx - Qx) - (Qx - Px) * (Ry - Qy),
  case Val == 0 of
    true ->
      0;
    false ->
      case Val > 0 of
        true ->
          1;
        _ ->
          2
      end
  end.

%% 確認 line(P1, Q1) 是否有跟 line(P2, Q2) 相交
-spec intersect(P1 :: point(), Q1 :: point(), P2 :: point(), Q2 :: point() ) -> boolean().
intersect(P1, Q1, P2, Q2) ->
  Orientation1 = orientation(P1, Q1, P2),
  Orientation2 = orientation(P1, Q1, Q2),
  Orientation3 = orientation(P2, Q2, P1),
  Orientation4 = orientation(P2, Q2, Q1),

  % general case
  case (Orientation1 /= Orientation2) and (Orientation3 /= Orientation4) of
    true ->
      true;
    _ ->
      %% Special Cases
      %% p1, q1 and p2 are colinear and p2 lies on segment p1q1
      case (Orientation1 == 0) and on_segment(P1, P2, Q1) of
        true ->
          true;
        _ ->
          %% p1, q1 and p2 are colinear and q2 lies on segment p1q1
          case (Orientation2 == 0) and on_segment(P1, Q2, Q1) of
            true ->
              true;
            _ ->
              %% p2, q2 and p1 are colinear and  p1 lies on segment p2q2
              case (Orientation3 == 0) and on_segment(P2, P1, Q2) of
                true ->
                  true;
                _ ->
                  %% p2, q2 and q1 are colinear and q1 lies on segment p2q2
                  case (Orientation4 == 0) and on_segment(P2, Q1, Q2) of
                    true ->
                      true;
                    _ ->
                      false
                  end
              end
          end
      end
  end.

-spec in_polygon(Polygon :: list(), P :: point() ) -> boolean().
in_polygon(Polygon, P) ->
  case length(Polygon) < 3 of
    true ->
      false;
    _ ->
      #point{x = _Px, y = Py} = P,
      %% 產生一個點,最後要跟 P 連成一條 Py 到 INF 的水平線段
      PInf = #point{x=?INF, y=Py},

      %% 計算 line(P, PInf) 跟所有多邊形的邊線的交點的數量
      %% Count intersections of the above line with sides of polygon

      % 分成第一個點, 跟其他點 兩個 list
      {PolygonHead, PolygonLast} = lists:split(1, Polygon),

%%      io:format("PolygonHead=~p, PolygonLast=~p~n", [PolygonHead, PolygonLast]),

      % 把 第一個點接到 PolygonLast 後面
      Polygon2 = lists:append(PolygonLast, PolygonHead),
      % 合併 Polygon, Polygon2 為新的 list, [{Polygon, Polygon2}]
      PolygonList = lists:zip(Polygon, Polygon2),

%%      io:format("PolygonList=~p~n", [PolygonList]),

      {CountRes, OnSegmentFlagRes, OnSegmentRes} = lists:foldl(fun({P1, P2}, {Count, OnSegmentFlag, OnSegment}) ->

%%        io:format("  lists:foldl P1=~p, P2=~p, intersect(P1, P2, P, PInf)=~p, orientation(P1, P, P2)=~p, on_segment(P1, P, P2)=~p~n", [P1, P2, intersect(P1, P2, P, PInf), orientation(P1, P, P2), on_segment(P1, P, P2)]),
        %% 判斷 (P1, P2), (P, PInf) 是否有交點
        case intersect(P1, P2, P, PInf) of
          true ->
            case orientation(P1, P, P2) == 0 of
              true ->
                %% 如果 P 跟 line(P1, P2) 共線,判斷是否 P 有在該線段上
                {Count, true, on_segment(P1, P, P2)};
              _ ->
                {Count + 1, OnSegmentFlag, OnSegment}
            end;
          _ ->
            {Count, OnSegmentFlag, OnSegment}
        end
                                                               end,
        {0, false, false}, PolygonList),

%%      io:format("CountRes=~p, OnSegmentFlagRes=~p, OnSegmentRes=~p, (CountRes rem 2)=~p~n", [CountRes, OnSegmentFlagRes, OnSegmentRes, CountRes rem 2]),
      %% 判斷交點數量是否為奇數
      %% Return true if count is odd, false otherwise
      case OnSegmentFlagRes of
        true ->
          OnSegmentRes;
        _ ->
          case (CountRes rem 2) == 1 of
            true ->
              true;
            _ ->
              false
          end
      end
  end.

test() ->
  P = #point{x = 1.0, y = 2.0},
  Q = #point{x = 2.0, y = 4.0},
  R = #point{x = 3.0, y = 6.0},
  S = #point{x = 4.0, y = 8.0},

  ResOnSegment = on_segment(P, Q, R),
  io:format("ResOnSegment=~p~n", [ResOnSegment]),

  ResOrientation = orientation(P, Q, R),
  io:format("ResOrientation=~p~n", [ResOrientation]),

  ResIntersect = intersect(P, Q, R, S),
  io:format("ResIntersect=~p~n", [ResIntersect]),
  ok.

test2() ->
  Polygon = [ #point{x = 0.0, y = 0.0}, #point{x = 10.0, y = 0.0}, #point{x = 10.0, y = 10.0}, #point{x = 0.0, y = 10.0} ],

  Res1 = in_polygon(Polygon, #point{x = 20.0, y = 20.0}),
  io:format("Res1=~p~n", [Res1]),

  Res2 = in_polygon(Polygon, #point{x = 5.0, y = 5.0}),
  io:format("Res2=~p~n", [Res2]),

  Res3 = in_polygon(Polygon, #point{x = -1.0, y = 10.0}),
  io:format("Res3=~p~n", [Res3]),

  %%%%%%%%%%%
  Polygon2 = [ #point{x = 0.0, y = 0.0}, #point{x = 5.0, y = 5.0}, #point{x = 5.0, y = 0.0} ],
  Res4 = in_polygon(Polygon2, #point{x = 3.0, y = 3.0}),
  io:format("Res4=~p~n", [Res4]),

  Res5 = in_polygon(Polygon2, #point{x = 5.0, y = 1.0}),
  io:format("Res5=~p~n", [Res5]),

  Res6 = in_polygon(Polygon2, #point{x = 8.0, y = 1.0}),
  io:format("Res6=~p~n", [Res6]),

  ok.

2020/08/10

詞向量 Word Embedding

文章本身是一種非結構化的資料,無法直接被計算。word representation 就是將這種訊息轉化為結構化的資訊,這樣就可以針對 word representation 計算,完成文章分類、情緒判斷等工作。

word representation 的方法很多,例如:

  1. one-hot representation

    例如: 貓、狗、牛、羊 用向量中一個欄位來表示

    貓:[1,0,0,0]
    狗:[0,1,0,0]
    牛:[0,0,1,0]
    羊:[0,0,0,1]

    缺點是沒有辦法表示出詞語之間的關係,另外因為向量中大部分都是 0,稀疏的向量,導致計算及儲存的效率都很低

  2. integer representation

    都以一個整數來表示每一個詞,將詞語的整數連接成 list,就是一句話

    貓:1
    狗:2
    牛:3
    羊:4

    缺點是沒有辦法表示出詞語之間的關係

  3. word embedding

    可用較低維度的向量表示詞語,不像one-hot 那麼長。詞意相近的詞,在向量空間中的距離比較接近

    有兩種主流的 word embedding 方法

    • word2vec

      2013 年由 google 的 Mikolov 提出,該演算法有兩種模式:利用前後文來預測目前的詞語,或是利用目前的詞語預測前後文

    • GloVe (Global Vector for Word Representation)

      延伸了 word2vec

word2vec

  • CBOW (Continuous Bag-of-Words Model)

    利用前後文來預測目前的詞語,相當於一句話中扣掉一個詞,猜這個詞是什麼。

  • Skip-gram (Continuous Skip-gram Model)

    利用目前的詞語預測前後文,相當於給一個詞,猜前面和後面可能出現什麼詞。

    

ref: Word2Vec 的兩種模型:CBOW 與 Skip-gram

優點:

  1. 通用性佳,適合用在多種 NLP 問題上
  2. 比舊的 word embedding 方法的向量維度小,計算速度比較快

缺點:

  1. 由於詞和向量是一對一的關係,所以無法處理多義詞的問題
  2. word2vec 是一種靜態的表示方式,通用性強,但無法針對特定任務做動態優化

window

以 「孔乙己 一到 店 所有 喝酒 的 人 便都 看著 他 笑」 這一句話為例,去掉停用字後,會得到

孔乙己 一到 店 所有 喝酒 人 看著 笑

以 「人」 這個單詞為例,window =1 時,就是該單詞前後 1 格的另一個單詞,這樣會得到這樣的結果

喝酒 人 看著

電腦就能知道「人」跟「喝酒」「看著」有關係。

window 用來定義 word2vec 文章分析時,單詞前後關係的距離。

gensim

gensim 是使用 google 釋出的 word2vec 模型的套件,可找到字的向量、相似字,計算向量之間的相似度,WMDistance 可計算兩個句子之間的相似度。

取得 wiki 文章資料

以下以 維基百科 wiki zh data 下載的 20200301 中文版資料 zhwiki-20200301-pages-articles.xml.bz2 1.8 GB 做測試,注意我們要的是以 pages-articles.xml.bz2 結尾的備份。

先安裝 gensim

pip3 install gensim

gensim 已經有提供了 WikiCorpus,可以快速取得 wiki 文章的標題及內容。執行以下程式,會產生一個 wiki_texts.txt 文字檔,裡面是所有 wiki_corpus.get_texts() 取得的文章內容。

# -*- coding: utf-8 -*-
## wiki_to_txt.py

import logging
import sys

from gensim.corpora import WikiCorpus

def main():

    if len(sys.argv) != 2:
        print("Usage: python3 " + sys.argv[0] + " wiki_data_path")
        exit()

    logging.basicConfig(format='%(asctime)s : %(levelname)s : %(message)s', level=logging.INFO)
    wiki_corpus = WikiCorpus(sys.argv[1], dictionary={})
    texts_num = 0

    with open("wiki_texts.txt",'w',encoding='utf-8') as output:
        for text in wiki_corpus.get_texts():
            output.write(' '.join(text) + '\n')
            texts_num += 1
            if texts_num % 10000 == 0:
                logging.info("已處理 %d 篇文章" % texts_num)

if __name__ == "__main__":
    main()

執行

python3 wiki_to_txt.py zhwiki-20200301-pages-articles.xml.bz2

執行結果

2020-03-30 15:27:59,410 : INFO : finished iterating over Wikipedia corpus of 356901 documents with 82295378 positions (total 3436353 articles, 97089149 positions before pruning articles shorter than 50 words)

斷詞

因為wiki 文章中,把簡體字跟繁體字混在一起了,先透過 OpenCC 進行簡體字轉繁體的處理

ref:

安裝 OpenCC

wget https://github.com/BYVoid/OpenCC/archive/ver.1.0.5.tar.gz -O opencc.1.0.5.tgz

tar -zxvf opencc.1.0.5.tgz
cd OpenCC-ver.1.0.5/

# 產生 Makefile
mkdir build
cd build

## CENTOS 執行以下命令
cmake -DCMAKE_INSTALL_PREFIX=/usr -DCMAKE_BUILD_TYPE=Release LE_GETTEXT:BOOL=ON  ..
## MAC 執行以下命令
# cmake -DCMAKE_INSTALL_PREFIX=/usr/local -DCMAKE_BUILD_TYPE=Release -D ENABLE_GETTEXT:BOOL=OFF  -DCMAKE_OSX_ARCHITECTURES=x86_64  ..

make
make install

sudo ln -s /usr/lib/libopencc.so.2 /usr/lib64/libopencc.so.2

## 測試
opencc --help
opencc --version

利用 opencc 將簡體字轉為繁體

opencc -i wiki_texts.txt -o wiki_zh_tw.txt -c s2tw.json

安裝 jieba

pip3 install jieba

測試

import jieba

seg_list = jieba.cut("我来到清华大学", cut_all=False)
print("Default Mode: " + "/ ".join(seg_list))

ithomeironman/day16NLP_Chinese/ 可下載一個繁體中文的字典 dict.txt.big,以及 停用字 stops.txt。將 jieba 改為使用繁體字典

# -*- coding: utf-8 -*-
## segment.py
import jieba
import logging

def main():

    logging.basicConfig(format='%(asctime)s : %(levelname)s : %(message)s', level=logging.INFO)

    # jieba custom setting.
    jieba.set_dictionary('jieba_dict/dict.txt.big')

    # load stopwords set
    stopword_set = set()
    with open('jieba_dict/stops.txt','r', encoding='utf-8') as stopwords:
        for stopword in stopwords:
            stopword_set.add(stopword.strip('\n'))

    output = open('wiki_seg.txt', 'w', encoding='utf-8')
    with open('wiki_zh_tw.txt', 'r', encoding='utf-8') as content :
        for texts_num, line in enumerate(content):
            line = line.strip('\n')
            words = jieba.cut(line, cut_all=False)
            for word in words:
                if word not in stopword_set:
                    output.write(word + ' ')
            output.write('\n')

            if (texts_num + 1) % 10000 == 0:
                logging.info("已完成前 %d 行的斷詞" % (texts_num + 1))
    output.close()

if __name__ == '__main__':
    main()

執行要花 30 分鐘

# python3 segment.py
2020-03-30 15:50:14,355 : DEBUG : Prefix dict has been built successfully.
......
2020-03-30 16:17:47,295 : INFO : 已完成前 350000 行的斷詞

訓練單詞向量

Word2Vec 有許多參數

gensim.models.word2vec.Word2Vec(sentences=None, size=100, alpha=0.025, window=5, min_count=5, max_vocab_size=None, sample=0.001, seed=1, workers=3, min_alpha=0.0001, sg=0, hs=0, negative=5, cbow_mean=1, hashfxn=<built-in function hash>, iter=5, null_word=0, trim_rule=None, sorted_vocab=1, batch_words=10000)

比較常用的是

  • sentences:這是要訓練的句子集合
  • size:這是訓練出的詞向量會有幾維
  • alpha:機器學習中的學習率,這東西會逐漸收斂到 min_alpha
  • sg:sg=1表示採用skip-gram,sg=0 表示採用cbow
  • window:能往左往右看幾個字的意思
  • workers:執行緒數目
  • min_count:若這個詞出現的次數小於min_count,那他就不會被視為訓練對象
# -*- coding: utf-8 -*-

import logging

from gensim.models import word2vec

def main():

    logging.basicConfig(format='%(asctime)s : %(levelname)s : %(message)s', level=logging.INFO)
    sentences = word2vec.LineSentence("wiki_seg.txt")
    model = word2vec.Word2Vec(sentences, size=250)

    #保存模型,供日後使用
    model.save("word2vec.model")

    #模型讀取方式
    # model = word2vec.Word2Vec.load("your_model_name")

if __name__ == "__main__":
    main()

模型測試

# -*- coding: utf-8 -*-

from gensim.models import word2vec
from gensim import models
import logging

def main():
    logging.basicConfig(format='%(asctime)s : %(levelname)s : %(message)s', level=logging.INFO)
    model = models.Word2Vec.load('word2vec.model')

    print("提供 3 種測試模式\n")
    print("輸入一個詞,則去尋找前一百個該詞的相似詞")
    print("輸入兩個詞,則去計算兩個詞的餘弦相似度")
    print("輸入三個詞,進行類比推理")

    while True:
        try:
            query = input()
            q_list = query.split()

            if len(q_list) == 1:
                print("相似詞前 100 排序")
                res = model.most_similar(q_list[0],topn = 100)
                for item in res:
                    print(item[0]+","+str(item[1]))

            elif len(q_list) == 2:
                print("計算 Cosine 相似度")
                res = model.similarity(q_list[0],q_list[1])
                print(res)
            else:
                print("%s之於%s,如%s之於" % (q_list[0],q_list[2],q_list[1]))
                res = model.most_similar([q_list[0],q_list[1]], [q_list[2]], topn= 100)
                for item in res:
                    print(item[0]+","+str(item[1]))
            print("----------------------------")
        except Exception as e:
            print(repr(e))

if __name__ == "__main__":
    main()

測試

籃球
相似詞前 100 排序
美式足球,0.6760541796684265
排球,0.6475502848625183
橄欖球,0.6430544257164001
男子籃球,0.6427032351493835
冰球,0.6138877272605896
棒球,0.6081532835960388
籃球隊,0.6004550457000732
足球,0.5992617607116699
.....

----------------------------
電腦 程式
計算 Cosine 相似度
0.5263175
----------------------------
衛生紙 啤酒
計算 Cosine 相似度
0.3263663
----------------------------
衛生紙 面紙
計算 Cosine 相似度
0.70471543
----------------------------

電腦 程式 電視
電腦之於電視,如程式之於
電腦系統,0.6098563075065613
程式碼,0.6063085198402405
軟體,0.5896543264389038
電腦程式,0.5740373730659485
終端機,0.5652117133140564
計算機程序,0.5597981810569763
除錯,0.554024875164032
計算機,0.549680769443512
作業系統,0.543748140335083
直譯器,0.5432565212249756
介面,0.5425338745117188
......

自然語言處理的應用

簡單

  • 拼寫檢查 ( Spell Checking )
  • 關鍵字搜索 ( Keyword Search )
  • 尋找同義詞 ( Finding Synonyms )

  • 從網頁和文檔解析信息 ( Parsing information from websites, documents, etc. )

複雜

  • 機器翻譯 ( Machine Translation )
  • 語義分析 ( Semantic Analysis )
  • 指代詞分析 ( Coreference ), 例如,”he” 和”it” 在文檔中指誰或什麼?
  • 問答系統 ( Question Answering )

References

詞向量詳解:從word2vec、glove、ELMo到BERT

Word2vec

詞嵌入 | Word embedding

自然語言處理入門- Word2vec小實作

讓電腦聽懂人話: 直觀理解 Word2Vec 模型

Gensim Word2Vec 簡易教學

產品標籤分群實作-Word2Vec

以 gensim 訓練中文詞向量

實作Tensorflow (5):Word2Vec

2020/08/03

中文斷詞

在中文自然語言處理NLP中,要對一堆文字詞語組成的文章進行分析, 分析前要先拆解文章,也就是斷詞,我們要分析的對象是詞語,而不是一個一個中文字,這跟英文完全不同,因為英文的斷詞就直接用標點符號、空白去區隔即可。

目前繁體中文斷詞系統有 中研院 CKIP 以及 jieba,在一些舊的文章中都提到 jieba 無法適當地處理繁體中文,而有替換繁體中文字典的改進作法,不過目前 jieba 已經到了 0.42 版,以下先了解官方的套件的功能,再看看需不需要修改繁體中文字典。

jieba 演算法

  • 基於前綴詞典實現高效的詞圖掃描,生成句子中漢字所有可能成詞情況所構成的有向無環圖 (DAG)
  • 採用了動態規劃查找最大概率路徑,找出基於詞頻的最大切分組合
  • 對於未登錄詞,採用了基於漢字成詞能力的 HMM 模型,使用了 Viterbi 算法

安裝

可直接用 pip 安裝,或是將 jieba source code 的 jieba 目錄放在目前的工作目錄,或是 site-packages 目錄中

如果要使用 paddle 的分詞語詞性標注功能,必須安裝 paddlepaddle-tiny

pip3 install paddlepaddle-tiny==1.6.1

先直接下載 source code 試試看

wget https://github.com/fxsjy/jieba/archive/v0.42.1.tar.gz -O jieba-0.41.1.tgz

tar zxvf jieba-0.41.1.tgz

virtual environemnt

virtualenv --system-site-packages /root/venv-jieba
source /root/venv-jieba/bin/activate

## 如果不使用 paddlepaddle,這兩個套件也可以不安裝
pip3 install numpy==1.16.4
pip3 install paddlepaddle-tiny==1.6.1

# 把 soruce code 中的 jieba 目錄移動到工作目錄中
mv ~/temp/download/jieba-0.41.1/jieba ~/temp

斷詞

有四種斷詞模式

  • 精確模式,試圖將句子最精確地切開,適合文本分析
  • 完整模式,把句子中所有的可以成詞的詞語都掃描出來, 速度非常快,但是不能解決歧義;
  • 搜索引擎模式,在精確模式的基礎上,對長詞再次切分,提高召回率,適合用於搜索引擎分詞。
  • paddle模式,利用PaddlePaddle深度學習框架,訓練序列標注(雙向GRU)網絡模型實現分詞。同時支持詞性標注。paddle模式使用需安裝 paddlepaddle-tiny

函式

  • jieba.cut 方法接受四個輸入參數: 需要分詞的字符串;cutall 參數用來控制是否採用全模式;HMM 參數用來控制是否使用 HMM 模型;usepaddle 參數用來控制是否使用paddle模式下的分詞模式,paddle模式採用延遲加載方式,通過enable_paddle接口安裝paddlepaddle-tiny,並且import相關代碼
  • jieba.cutforsearch 方法接受兩個參數:需要分詞的字符串;是否使用 HMM 模型。該方法適合用於搜索引擎構建倒排索引的分詞,粒度比較細
  • 待分詞的字符串可以是 unicode 或 UTF-8 字符串、GBK 字符串。注意:不建議直接輸入 GBK 字符串,可能無法預料地錯誤解碼成 UTF-8
  • jieba.cut 以及 jieba.cutforsearch 返回的結構都是一個可迭代的 generator,可以使用 for 循環來獲得分詞後得到的每一個詞語(unicode),或者用 jieba.lcut 以及 jieba.lcutforsearch 直接返回 list
  • jieba.Tokenizer(dictionary=DEFAULT_DICT) 新建自定義分詞器,可用於同時使用不同詞典。jieba.dt 為默認分詞器,所有全局分詞相關函數都是該分詞器的映射。

# encoding=utf-8
import jieba

print()
print("完整模式:")
seg_list = jieba.cut("肺炎疫情的挑戰日益嚴峻,新竹清華大學自農曆年起陸續已採取了量測體溫等全面的防疫措施。", cut_all=True)
print("Full Mode: " + "/ ".join(seg_list))  # 全模式

print()
print("精確模式:")
seg_list = jieba.cut("肺炎疫情的挑戰日益嚴峻,新竹清華大學自農曆年起陸續已採取了量測體溫等全面的防疫措施。", cut_all=False)
print("Default Mode: " + "/ ".join(seg_list))  # 精確模式

print()
print("預設是精確模式:")
seg_list = jieba.cut("肺炎疫情的挑戰日益嚴峻,新竹清華大學自農曆年起陸續已採取了量測體溫等全面的防疫措施。")  # 預設是精確模式
print(", ".join(seg_list))

print()
print("搜索引擎模式:")
seg_list = jieba.cut_for_search("肺炎疫情的挑戰日益嚴峻,新竹清華大學自農曆年起陸續已採取了量測體溫等全面的防疫措施。")  #
print(", ".join(seg_list))

print()
print("Paddle Mode:")
jieba.enable_paddle()# 啓動paddle模式。 0.40版之後開始支持,早期版本不支持
strs=["肺炎疫情的挑戰日益嚴峻,新竹清華大學自農曆年起陸續已採取了量測體溫等全面的防疫措施。","乒乓球拍賣完了","新竹清華大學"]
for str in strs:
    seg_list = jieba.cut(str,use_paddle=True) # 使用paddle模式
    print("Paddle Mode: " + '/'.join(list(seg_list)))

執行結果

完整模式:
Building prefix dict from the default dictionary ...
Loading model from cache /tmp/jieba.cache
Loading model cost 0.433 seconds.
Prefix dict has been built successfully.
Full Mode: 肺炎/ 疫情/ 的/ 挑/ 戰/ 日益/ 嚴/ 峻/ ,/ 新竹/ 清/ 華/ 大/ 學/ 自/ 農/ 曆/ 年/ 起/ 陸/ 續/ 已/ 採/ 取/ 了/ 量/ 測/ 體/ 溫/ 等/ 全面/ 的/ 防疫/ 措施/ 。

精確模式:
Default Mode: 肺炎/ 疫情/ 的/ 挑戰/ 日益/ 嚴峻/ ,/ 新竹/ 清華大學/ 自農/ 曆/ 年/ 起/ 陸續/ 已/ 採取/ 了/ 量/ 測體/ 溫/ 等/ 全面/ 的/ 防疫/ 措施/ 。

預設是精確模式:
肺炎, 疫情, 的, 挑戰, 日益, 嚴峻, ,, 新竹, 清華大學, 自農, 曆, 年, 起, 陸續, 已, 採取, 了, 量, 測體, 溫, 等, 全面, 的, 防疫, 措施, 。

搜索引擎模式:
肺炎, 疫情, 的, 挑戰, 日益, 嚴峻, ,, 新竹, 清華大學, 自農, 曆, 年, 起, 陸續, 已, 採取, 了, 量, 測體, 溫, 等, 全面, 的, 防疫, 措施, 。

Paddle Mode:
W0327 11:45:31.179752 21561 init.cc:157] AVX is available, Please re-compile on local machine
Paddle enabled successfully......
Paddle Mode: 肺炎/疫情/的/挑戰/日益/嚴/峻,新竹清華大學自農曆年起陸續/已/採取/了/量測體溫/等/全面/的/防疫/措施/。
Paddle Mode: 乒乓/球拍/賣/完/了
Paddle Mode: 新竹/清華大學
  • jieba.lcut() lcut(),意思跟cut()是一樣的,只是返回的型態變成list

# encoding=utf-8
import jieba

print()
print("完整模式:")
seg_list = jieba.lcut("肺炎疫情的挑戰日益嚴峻,新竹清華大學自農曆年起陸續已採取了量測體溫等全面的防疫措施。", cut_all=True)
print("Full Mode: ", seg_list)  # 全模式


print()
print("精確模式:")
seg_list = jieba.lcut("肺炎疫情的挑戰日益嚴峻,新竹清華大學自農曆年起陸續已採取了量測體溫等全面的防疫措施。", cut_all=False)
print("Default Mode: ", seg_list)  # 精確模式

print()
print("預設是精確模式:")
seg_list = jieba.lcut("肺炎疫情的挑戰日益嚴峻,新竹清華大學自農曆年起陸續已採取了量測體溫等全面的防疫措施。")  # 預設是精確模式
print(seg_list)

print()
print("搜索引擎模式:")
seg_list = jieba.lcut_for_search("肺炎疫情的挑戰日益嚴峻,新竹清華大學自農曆年起陸續已採取了量測體溫等全面的防疫措施。")  #
print(seg_list)

print()
print("Paddle Mode:")
jieba.enable_paddle()# 啓動paddle模式。 0.40版之後開始支持,早期版本不支持
strs=["肺炎疫情的挑戰日益嚴峻,新竹清華大學自農曆年起陸續已採取了量測體溫等全面的防疫措施。","乒乓球拍賣完了","新竹清華大學"]
for str in strs:
    seg_list = jieba.lcut(str,use_paddle=True) # 使用paddle模式
    print("Paddle Mode: ", seg_list)

自訂詞典

  • 雖然 jieba 有新詞識別能力,但是自行添加新詞可以保證更高的正確率

  • 用法: jieba.loaduserdict(filename)

    file_name 為文件類對象或自定義詞典的路徑

  • 詞典格式和 dict.txt 一樣,一個詞佔一行;每一行分三部分:詞語、詞頻(可省略)、詞性(可省略),用空格隔開,順序不可顛倒。file_name 若為路徑或二進制方式打開的文件,則文件必須為 UTF-8 編碼。

  • 詞頻省略時使用自動計算的方式處理,能保證分出該詞的詞頻。

  • 使用 addword(word, freq=None, tag=None) 和 delword(word) 可在程序中動態修改詞典。

  • 使用 suggest_freq(segment, tune=True) 可調節單個詞語的詞頻,使其能(或不能)被分出來。

  • 注意:自動計算的詞頻在使用 HMM 新詞發現功能時可能無效。

在專案路徑下新增一個檔案叫做:userdict.txt

內容如下:

農曆年
量測
體溫
日益嚴峻

可在程式一開始,就載入自訂詞典

jieba.load_userdict('userdict.txt')

執行結果

精確模式:
Default Mode:  ['肺炎', '疫情', '的', '挑戰', '日益嚴峻', ',', '新竹', '清華大學', '自', '農曆年', '起陸續', '已', '採取', '了', '量測', '體溫', '等', '全面', '的', '防疫', '措施', '。']

可動態調整詞語的頻率

jieba.suggest_freq(('陸續'), True)
print()
print("精確模式2:")
seg_list = jieba.lcut("肺炎疫情的挑戰日益嚴峻,新竹清華大學自農曆年起陸續已採取了量測體溫等全面的防疫措施。", cut_all=False)
print("Default Mode: ", seg_list)  # 精確模式

執行結果

精確模式2:
Default Mode:  ['肺炎', '疫情', '的', '挑戰', '日益嚴峻', ',', '新竹', '清華大學', '自', '農曆年', '起', '陸續', '已', '採取', '了', '量測', '體溫', '等', '全面', '的', '防疫', '措施', '。']

也可以直接替換詞典

ithomeironman/day16NLP_Chinese/ 可下載一個繁體中文的字典 dict.txt.big

# encoding=utf-8
import jieba

jieba.set_dictionary('dict.txt.big')
# with open('stops.txt', 'r', encoding='utf8') as f:
#     stops = f.read().split('\n')

print()
print("精確模式:")
seg_list = jieba.lcut("肺炎疫情的挑戰日益嚴峻,新竹清華大學自農曆年起陸續已採取了量測體溫等全面的防疫措施。", cut_all=False)
print("Default Mode: ", seg_list)

關鍵詞抽取

TF-IDF 方法

# jieba.analyse.extract_tags(sentence, topK=20, withWeight=False, allowPOS=())
# sentence 為待提取的文本
# topK 為返回幾個 TF/IDF 權重最大的關鍵詞,默認值為 20
# withWeight 為是否一並返回關鍵詞權重值,默認值為 False
# allowPOS 僅包括指定詞性的詞,默認值為空,即不篩選
# jieba.analyse.TFIDF(idf_path=None) 新建 TFIDF 實例,idf_path 為 IDF 頻率文件

TextRank 方法

jieba.analyse.textrank(sentence, topK=20, withWeight=False, allowPOS=('ns', 'n', 'vn', 'v')) 直接使用,接口相同,注意默認過濾詞性。
jieba.analyse.TextRank() 新建自定義 TextRank 實例
# -*- coding: utf-8 -*-
import jieba
import jieba.analyse

jieba.set_dictionary('dict.txt.big')

print()
text = '肺炎疫情的挑戰日益嚴峻,新竹清華大學自農曆年起陸續已採取了量測體溫等全面的防疫措施。'
tags = jieba.analyse.extract_tags(text, topK=10)
print(tags)
# ['挑戰', '嚴峻', '清華大學', '農曆年', '陸續', '採取', '測體溫', '防疫', '新竹', '肺炎']

print()
print('textrank:')
for x, w in jieba.analyse.textrank(text, withWeight=True):
    print('%s %s' % (x, w))

# textrank:
# 日益 1.0
# 全面 1.0
# 肺炎 0.6631715416020616
# 防疫 0.6631715416020616
# 疫情 0.6605033585768562
# 措施 0.6605033585768562
# 新竹 0.3607120276929184
# 了量 0.3607120276929184

詞性標注

ref: 彙整中文與英文的詞性標註代號:結巴斷詞器與FastTag / Identify the Part of Speech in Chinese and English

結巴預設會將標點符號標示為「x」,而不是「w」。而且英文會被標示為「eng」

# -*- coding: utf-8 -*-
import jieba
import jieba.posseg as pseg

text = '肺炎疫情的挑戰日益嚴峻,新竹清華大學自農曆年起陸續已採取了量測體溫等全面的防疫措施。'
seg_list = pseg.lcut(text)
# print("Default Mode: ", seg_list)

for word, flag in seg_list:
    print("", word, " : ", flag)

執行結果

 肺炎  :  n
 疫情  :  n
 的  :  uj
 挑戰  :  vn
 日益  :  n
 嚴峻  :  a
 ,  :  x
 新竹  :  ns
 清華大學  :  nt
 自  :  p
 農  :  ng
 曆  :  zg
 年  :  q
 起  :  v
 陸  :  nr
 續  :  v
 已  :  d
 採  :  v
 取  :  v
 了  :  ul
 量  :  n
 測  :  v
 體  :  ng
 溫  :  v
 等  :  u
 全面  :  n
 的  :  uj
 防疫  :  vn
 措施  :  n
 。  :  x

word cloud

參考這篇文章中文自然語言處理基礎 以及資源

下載檔案:

cloud_mask7.png
sumsun.ttf
stops.txt

安裝其它套件

pip3 install collections
pip3 install wordcloud
pip3 install matplotlib
yum -y install python-imaging
import jieba
jieba.set_dictionary('dict.txt.big')  # 如果是使用繁體文字,請記得去下載繁體字典來使用
with open('stops.txt', 'r', encoding='utf8') as f:
    stops = f.read().split('\n')

text = "肺炎疫情的挑戰日益嚴峻,新竹清華大學自農曆年起陸續已採取了量測體溫等全面的防疫措施。"

from collections import Counter
from wordcloud import WordCloud
from matplotlib import pyplot as plt

stops.append('\n')  ## 換行符號,加入停用字中,可以把它拿掉
stops.append('\n\n')
terms = [t for t in jieba.cut(text, cut_all=True) if t not in stops]

sorted(Counter(terms).items(), key=lambda x:x[1], reverse=True)  ## 這個寫法很常出現在Counter中,他可以排序,list每個item出現的次數。


plt.clf()
wordcloud = WordCloud(font_path="simsun.ttf")  ##做中文時務必加上字形檔
wordcloud.generate_from_frequencies(frequencies=Counter(terms))
plt.figure(figsize=(15,15))
plt.imshow(wordcloud, interpolation="bilinear")
plt.axis("off")
plt.savefig('cloud1.png')


from PIL import Image
import numpy as np

alice_mask = np.array(Image.open("cloud_mask7.png"))  ## 請更改cloud_mask7.png路徑
wc = WordCloud(background_color="white", max_words=2000, mask=alice_mask, font_path="simsun.ttf")
wc.generate_from_frequencies(Counter(terms))  ## 請更改Counter(terms)

wc.to_file("cloud2.png")  ##如果要存檔,可以使用

# plt.clf()
# plt.imshow(wc, interpolation='bilinear')
# plt.axis("off")
# plt.figure()
# plt.imshow(alice_mask, cmap=plt.cm.gray, interpolation='bilinear')
# plt.axis("off")
# plt.savefig('cloud2.png')

中文斷詞.png

References

pypi jieba

jieba 原始 github

python-11-利用jieba實現中文斷詞

NLP 中文斷詞最方便的開源工具之一 —— Jieba

Python-知名Jieba中文斷詞工具教學

2017

結巴中文斷詞台灣繁體版本 APCLab

結巴中文斷詞台灣繁體版本 ldkrsi

2020/07/20

馬可夫鏈 Markov Chain

有許多自然或社會現象,其發展或演變會隨時間進行而呈現出幾種可能的狀態,稱為狀態空間,這種隨機變化的過程稱為隨機過程 stochastic process。如果記錄狀態的時間點是離散的,就稱為離散時間隨機過程, Markov Process 或是高速公路每天的意外事件數量就是一種離散時間隨機過程,如果記錄狀態的時間點是連續的,就稱為連續時間隨機過程,布朗運動 Brownian Motion 或是在銀行等待服務的人數就是一種連續時間隨機過程。

如果時間與狀態都是離散的 Markov Process 就稱為 Markov Chain。Markov Chain 是隨機變數 \(X_0, X_1, X_2, X_3, ...\) ,也就是 \({X_n, n>=0}\) 的一個數列,其中 \(X_0, X_1, ... X_{n-1}\) 都是過去時間的狀態, \(X_n\) 是目前的狀態, \(X_{n+1}, X_{n+2}, ...\) 是未來的狀態,每一個狀態可轉移到其他狀態,轉移的機率總和為 1,因此目前的狀態都是由前一個狀態以及轉移機率來決定的。

討論 Markov Chain 有兩個先決條件

  1. 目前的狀態都是由前一個狀態的轉移機率來決定的
  2. 在任何時間點,系統的事件只存在於某一種狀態中

Makov Chain 的表示方式

  1. transition diagram

  2. transition matrix 轉移矩陣

    \( \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0.5 & 0.5 \\ 0.6 & 0 & 0.4 \end{bmatrix} \)

  3. tree diagram

transition matrix

為方便計算,以 transition matrix 作為表示方式

\( T = \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} \)

選定某一個狀態為起點,根據轉移矩陣改變狀態,當狀態不斷地改變,會產生一個狀態變化的路徑。如考慮所有的狀況,可計算出每一個路徑的發生機率,這些路徑就稱為 Markov Chain。


問題:如果一開始的狀態在 \(a_1\) ,三天後,在 \(a_2\) 的機率為?

一天後: \( \begin{bmatrix} 1 & 0 & 0 \end{bmatrix}\begin{bmatrix} 0 & 1 & 0 \\ 0 & 0.5 & 0.5 \\ 0.6 & 0 & 0.4 \end{bmatrix} = \begin{bmatrix} 0 & 1 & 0 \end{bmatrix} \)

兩天後: \( \begin{bmatrix} 0 & 1 & 0 \end{bmatrix}\begin{bmatrix} 0 & 1 & 0 \\ 0 & 0.5 & 0.5 \\ 0.6 & 0 & 0.4 \end{bmatrix} = \begin{bmatrix} 0 & 0.5 & 0.5 \end{bmatrix} \)

三天後: \( \begin{bmatrix} 0 & 0.5 & 0.5 \end{bmatrix}\begin{bmatrix} 0 & 1 & 0 \\ 0 & 0.5 & 0.5 \\ 0.6 & 0 & 0.4 \end{bmatrix} = \begin{bmatrix} 0.3 & 0.25 & 0.45 \end{bmatrix} \)

如果直接一次計算:\( \begin{bmatrix} 1 & 0 & 0 \end{bmatrix} T^3 = \begin{bmatrix} 1 & 0 & 0 \end{bmatrix}(\begin{bmatrix} 0 & 1 & 0 \\ 0 & 0.5 & 0.5 \\ 0.6 & 0 & 0.4 \end{bmatrix})^3 = \begin{bmatrix} 0.3 & 0.25 & 0.45 \end{bmatrix} \)

結果為 0.25

吸收馬可夫鏈

在 Markov Process 中,當事件進入某個狀態時,就一直停留在該狀態,稱該狀態為吸收狀態 absorbing state

如果 Markov Chain 有下列兩種特性,就稱為 Absorbing Markov Chain

  1. 至少存在一個吸收狀態
  2. 事件在任何非吸收狀態時,經過有限時間後,就會進入 absorbing state

ex:

\( \begin{bmatrix} 0.2 & 0.3 & 0 & 0.5 \\ 0 & 0.4 & 0 & 0.6 \\ 0.5 & 0 & 0.5 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \) 是一種 Absorbing Markov Chain

\( \begin{bmatrix} 0.2 & 0 & 0.4 & 0.4 \\ 0 & 1 & 0 & 0 \\ 0.1 & 0 & 0.5 & 0.4 \\ 0.5 & 0 & 0.3 & 0.2 \end{bmatrix} \) 不是一種 Absorbing Markov Chain,雖然 \(a_2\) 是 absorbing state,但是其他三種狀態,都不能到達 \(a_2\),因此這不是 Absorbing Markov Chain

醫院是一種 absorbing markov chain 的實例,如果將進入醫院住院的病患,分為住院可走動,住院臥床,痊癒出院,死亡四種狀態,這四種狀態可以構成一個 absorbing markov chain,其中痊癒出院即死亡,都分別為 absorbing state

工廠品管商品,生產線製造的新品(狀態1),通過品管後為合格品(狀態2),失敗的商品會進入再造品(狀態3),再造品可再重新加工,變成修復品(狀態4),重新加工失敗最終成為報廢品(狀態5)

\( \begin{matrix} \begin{matrix} \quad\quad\quad\quad & 合格品 & 報廢品 & 新品 & 再造品 & 修復品 \end{matrix} \\ \begin{matrix} 合格品 \\ 報廢品 \\ 新品 \\ 再造品 \\ 修復品 \end{matrix} \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0.7 & 0 & 0 & 0.3 & 0 \\ 0.9 & 0 & 0 & 0 & 0.1 \\ 0.6 & 0.4 & 0 & 0 & 0 \end{bmatrix} \end{matrix} \)

m 階馬可夫鏈

前面討論的 Markov Chain,用前一個單位時間的狀態可以預測下一個狀態,這稱為一階馬可夫鏈 first-order Markov Chain,如果需要前 m 個狀態,才能預測下一個狀態,就稱為 m 階馬可夫鏈。當然也有零階馬可夫鏈,也就是一般的機率分佈模型。

通常討論的馬可夫鏈都是一階馬可夫鏈,最重要的是無記憶性,因為系統不需要記錄前面經過的狀態,只需要目前的狀態,就可以預測下一個狀態。

HMM: Hidden Markov Model

當系統存在一些外部看不見的內部狀態,隨著系統內部狀態不同,表現出來的外部狀態也不同,我們可以直接觀察到外部狀態。

以買飲料為例,假設目前有三種飲料可以選擇:可樂、茶、牛奶,我們只能觀察到使用者每天早上9:00吃早餐,正在喝哪一種飲料,但是我們不知道他今天是在 7-11、FamilyMart 或是 HiLife 買到的。便利商店是內部狀態,三種飲料是外部狀態。

內部狀態的變化可用一個 transition matrix 表示

\( A= \begin{bmatrix} 0.7 & 0.1 & 0.2 \\ 0.3 & 0.2 & 0.5 \\ 0.4 & 0.3 & 0.3 \end{bmatrix} \)

不同的內部狀態,會產生不同的外部狀態,也可以用另一個 matrix 表示

\( B = \begin{bmatrix} 0.2 & 0.2 & 0.6 \\ 0.3 & 0.3 & 0.4 \\ 0.5 & 0.3 & 0.2 \end{bmatrix} \)

如果以 h 為內部狀態, \(h^{t}\) 為 t 單位時間的內部狀態, y 為外部狀態,\(y^{t}\) 是 t 單位時間的外部狀態,則

\(h^{t} = h^{t-1}A\)

\(y^{t} = h^{t}B\)

在實務上的問題,我們只知道外部狀態的觀察值,內部狀態被隱藏無法得知,也就是我們不會知道 transition matrix 實際的機率值,也就是我們不知道 A 與 B 的數值,只有一堆數據資料。

例如我們目前有甲乙丙三個人半年來每一週買飲料的歷史資料,當給定某一週買飲料的序列清單時,我們可以根據歷史資料,找出每一個人的購買習慣,也就是計算出 transition matrix,訓練模型,最後利用這個模型,將新的那一週的資料,分類判斷可能是哪一個人購買的。

在已知觀察序列,但未知轉移矩陣的條件下,找出其對應的的所有可能狀態序列,以及其中最大機率之狀態序列。

HMM 是有關時間序列的數據資料建立的模型,也是所有時間模型(ex: RNN, NLP, RL)的基礎。

HMM 三個基本問題

定義HMM參數為 λ,S 為狀態值集合, O 為觀察值集合,初始機率 π (Initial Probability)、轉移機率 A (Transition Probability)、發射機率 B (Emission Probability),因此 λ = ( π, A, B )

  1. Evaluation Problem

    評估問題最簡單,將所有機率做乘積就可以得到觀測序列發生的機率。

    根據觀察序列O及已知的模型參數λ,得出發生此序列的機率 P(O|λ)

    可用 Forward-backward Algorithm 降低計算複雜度

  2. Decoding Problem

    根據觀察序列O及模型參數λ,得出最有可能的隱藏狀態序列(hidden state sequence),也就是要求得哪一種狀態序列發生的機率最大,求給定觀測序列條件概率P(I|Q,λ)最大的狀態序列I

    可用 Viterbi Algorithm

    Viterbi算法是求最短路徑的過程,由於不知道每個觀測狀態的隱藏狀態是什麼,所以初始化時,每種隱藏狀態都有可能,計算在每種隱藏狀態下的整個觀測序列發生的機率,機率最大的那條路徑就是每個觀測狀態對應的唯一隱藏狀態。

  3. Learning Problem

    學習問題最複雜

    在已知隱藏狀態與觀察狀態的條件下,根據觀察序列O,找到最佳的參數 狀態轉移矩陣 A、觀察轉移矩陣 B、起始值機率 π,調整模型參數 λ = [A, B, π] 使得 P(O|λ) 最大化

    可用 EM Algorithm

語音辨識 HMM

語音辨識 HMM

語音辨識的目標,是事先使用大量語料,對每一個辨識語句建立一個 DHMM,然後當我們取得使用者輸入的測試語句後,即對測試語句進行端點偵測及 MFCC 特徵擷取,然後再將 MFCC 送至各個辨識語句所形成的 DHMM,算出每一個模型的機率值,機率值最大的 DHMM,所對應的辨識語句就是此系統的辨識結果。每一個音框的特徵向量都是39維的 MFCC

建立 0~9 數字語音辨識的方法

  1. 建立 0, 1, 2 ~ 9 的 HMMs
  2. 利用 Viterbi Decoding 找到新語句中,計算針對每一個 HMM 發音音素順序的序列的發生機率,也就是計算輸入語句和每個模型的相似度
  3. 最大機率值的 HMM,就是預測的 digit

進行 DHMM 的設計時,我們必須對每一個數字建立一個 HMM 模型,以數字「九」的發音為例,相關的 HMM 模型可以圖示如下:

把「九」的發音切分成三個「狀態」(States),分別代表 ㄐ、ㄧ、ㄡ 的發音,每一個狀態可以包含好幾個音框,而每一個音框隸屬於一個狀態的程度,是用一個機率值來表示,稱為「狀態機率」(State Probability)。此外,當一個新的音框進來時,我們可以將此音框留在目前的狀態,也可以將此音框送到下一個狀態,這些行為模式也都是用機率來表示,統稱為「轉移機率」(Transition Probability),其中我們使用「自轉移機率」(Self-transition Probability)來代表新音框留在目前狀態的機率,而用「次轉移機率」(Next-transition Probability)來代表新的音框跳到下一個狀態的機率。

由於每一個音框都會轉成一個語音特徵向量,我們首先使用 VQ 來將這些特徵向量轉成符號(Symbols),換句話說,每一個音框所對應的 VQ 群聚的索引值(Index),即是此音框的符號。若以數學來表示,k = O(i) 即代表 frame i 被分到 cluster k。

定義一些常用的數量:

  • frameNum:一段語音所切出的音框個數。
  • dim:每一個音框所對應的語音特徵向量的維度(例如基本的 MFCC 有 12 維)。
  • stateNum:一個 DHMM 的狀態個數
  • symbolNum:符號的個數,也就是 VQ 後的群數

HMM 的參數,就是指「狀態機率」和「轉移機率 Transition Probability」,說明如下:

  • 我們通常以矩陣 A 來表示轉移機率,其維度是 stateNum x stateNum,其中 A(i, j) 即是指由狀態 i 跳到狀態 j 的機率值。例如在上圖中,由狀態 1 跳到狀態 2 的機率是 0.3,因此 A(1, 2) = 0.3。一般而言,A(i, j) 滿足下列條件:
    • 在進行轉移時,我們只允許搜尋路徑跳到下一個相鄰的狀態,因此 A(i, j) = 0 if j≠i and j≠i+1。
    • 對於某一個狀態而言,所有的轉移機率總和為 1,因此 A(i, i) + A(i, i+1) = 1。
  • 我們通常以矩陣 B 來表示狀態機率,其維度是 symbolNum x stateNum,B(k,j) 即是指符號 k 隸屬於狀態 j 的機率值。換句話說,B 定義了由「符號」到「狀態」的機率,因此給定第 i 個音框,我們必須先由 O(i) 找出這個音框所對應的符號 k,然後再由 B(k,j) 找出此音框屬於狀態 j 的機率。

假設我們已經知道「九」的 HMM 所包含相關的參數 A 和 B,此時當我們錄音得到一段語音,我們如何算出來此語音隸屬於「九」的機率或程度?如何計算輸入語句和每個模型的相似度?換句話說,我們如何將每個音框分配到各個狀態之中,使得我們能夠得到整段語音最高的機率值?最常用的方法稱為 Viterbi Decoding,這是基於「動態規劃」(Dynamic Programming)的方法,可用數學符號定義如下:

  1. 目標函數:定義 D(i, j) 是 t(1:i) 和 r(1:j) 之間的最大的機率,t(1:i) 是前 i 個音框的特徵向量所成的矩陣,r(1:j) 則是由前 j 個狀態所形成的 DHMM,對應的最佳路徑是由 (1, 1) 走到 (i, j)。
  2. 遞迴關係:D(i, j) = B(O(i), j) + max{D(i-1, j)+A(j, j), D(i-1, j-1)+A(j-1, j)}
  3. 端點條件:D(1, j) = p(1, j) + B(O(1), j), j = 1 ~ stateNum
  4. 最後答案:D(m, n)

在上述數學式中,我們所用的所有機率都是「對數機率」(Log Probabilities),所以原來必須連乘的數學方程式,都變成了連加,具有下列好處:

  • 以加法來取代乘法,降低計算的複雜度。
  • 避開了由於連乘所造成的數值誤差。

假設轉移機率 A 和狀態機率 B 已知,那麼經由 Viterbi Decoding,我們可以找到最佳路徑(即是最佳分配方式),使整段路徑的機率值為最大,此機率值及代表輸入語音隸屬於此 DHMM 的機率,若是機率越高,代表此輸入語音越可能是「九」。

有關於 B(O(i),j)的定義,都是「音框 i 隸屬於狀態 j 的機率」,但是在 DHMM 和 CHMM 的算法差異很大,說明如下:

  • 在 DHMM,若要計算 B(O(i),j),我們必須先找到音框 i 所對應的符號 O(i),然後在經由查表法查到此符號 O(i) 屬於狀態 j 的機率是 B(O(i), j)。
  • 在 CHMM,B(O(i),j) 是由一個連續的機率密度函數所定義。

如何經由大量語料來估測每個 HMM 模型的最佳參數值 A 和 B?

首先,我們必須先定義什麼是「最佳參數」:對某一個特定模型而言,最佳參數值 A 和 B 應能使語料產生最大的對數機率總和。這個定義和最大似然估計(maximum likelihood estimation,縮寫為MLE) 完全相同,也非常合理。

計算最佳參數的方法,稱為 Re-estimation,其原理非常類似 batch k-means (Forgy's method) 的方法,先猜 A 和 B 的值,再反覆進行 Viterbi Decoding,然後再重新估算 A 和 B 的值,如此反覆計算,我們可以證明在疊代過程中,機率總和會一路遞增,直到逼近最佳值。(但是,就如同 k-means Clustering,我們並無法證明所得到的機率值是 Global Maximum,因此在訓練的過程中,可以使用不同的啟始參數,看看是否在反覆疊代後,能夠得到更好的機率總和。)求取參數的方法說明如下:

  1. 將所有的訓練語句切出音框,並將每一個音框轉成語音特徵向量,例如 39 維的 MFCC。

  2. 對所有語音特徵向量進行向量量化編碼(Vector Quantization, VQ),找出每一個向量所對應的 symbol(此即為對應的中心點或codeword)

  3. 先猜一組 A 和 B 的啟始值。如果沒有人工的標示資料,我們可以使用簡單的「均分法」,示意圖如下:

  4. 反覆進行下列兩步驟,直到收斂

    • Viterbi decoding:在 A 和 B 固定的情況下,利用 Viterbi decoding,找出 n 句「九」的語料所對應的 n 條最佳路徑。
    • Re-estimation:利用這 n 條最佳路徑,重新估算 A 和 B。

估算 A 的方法

\( A= \begin{bmatrix} 3/4 & 1/4 & 0 \\ 0 & 4/5 & 1/5 \\ 0 & 0 & 1 \end{bmatrix} \)

估算 B 的方法:

\( B= \begin{bmatrix} 1/4 & 3/5 & 0 \\ 1/4 & 2/5 & 2/6 \\ 2/4 & 0 & 0 \\ 0 & 0 & 4/6 \end{bmatrix} \)

假設我們的目標函數可以寫成 P(A, B, Path),則上述求參數的方法會讓 P(A, B, Path)逐次遞增,直到收斂,因為:

  1. 在 A, B 固定時,Viterbi decoding 會找出最佳的路徑,使 P(A, B, Path) 有最大值
  2. 在路徑(Path)固定時,Re-estimation 會找出最佳的 A, B 值以使 P(A, B, Path)有最大值

References

隨機過程.pdf

Markov Chain

Markov Model

Hidden Markov Model (part 1)

Markov chain 及 HMM

HMM(Hidden Markov Model)簡介

機器不學習:HMM模型解析

Hidden Markov Model (part 3)

hmmDiscrete_chinese

即時語音辨識系統