import numpy as np
import matplotlib.pyplot as plt
from copy import deepcopy

# Add Gaussian noise
def add_noise(x, mean=0, deviation=R, seed=None):

    if seed is not None:
        np.random.seed(seed=seed)

    shape = x.shape
    noise = np.random.normal(loc=mean, scale=deviation, size=shape)
    x_ = x + noise

    return x_

# Measurement time interval
dt = 1.

# Actual positions and velocity

v_ = 10. #deviation given to add_noise function in setting true_velocity

true_velocity = add_noise(np.ones((100,1),dtype=float), deviation=v_,seed=0)
true_positions = np.cumsum(true_velocity * dt)
true_positions = np.reshape(true_positions, [100,1])

true_x = np.concatenate((true_positions, true_velocity),axis=1)

#Observation

R = 1. #covariance of measurement noise
observed_x = add_noise(true_positions, deviation=R, seed=0)

#Initialization

z = deepcopy(observed_x)
x_init = true_x[0]
P_init = np.array([[0., 0.], [0., 0.]])
A = np.array([[1., dt], [0., 1.]])
Q = [0.1, 0.1] #process noise covariance
H = np.array([1., 0.]) #Observation Matrix: Extract observed values

#Kalman Filter
def kalman_filter(z, x_init, P_init, Q, R):
    
    #init_step
    x_filtered = np.array([x_init])
    P_filtered = np.array([P_init])   
    I = np.identity((len(x_init)))

    for t in range(1, len(z)):
        #Time Update(Predict)
        x = np.dot(A, x_filtered[t-1])
        P = np.dot(np.dot(A, P_filtered[t-1]), P_filtered[t-1].T) + Q

        #Measurement Update(Correct)
        K = np.dot(P, H.T) / (np.dot(np.dot(H, P), H.T) + R)
        y = z[t] - np.dot(H, x)
        x = x + (K * y)
        P = np.dot((I - np.dot(K, H)), P)

        #append filtered values
        x_filtered = np.append(x_filtered, np.array([x]), axis=0)
        P_filtered = np.append(P_filtered, np.array([P]), axis=0)
    
    return x_filtered, P_filtered

filtered_x = kalman_filter(z, x_init, P_init, Q=Q, R=R)[0]

plt.plot(true_x[:,0], 'r--', label='Actual Positions')
plt.plot(true_x[:,1], 'b--', label='Actual Velocity')
plt.plot(observed_x[:,0], 'y', label='Measured Positions')
plt.plot(filtered_x[:,0], 'gray', label='Estimated Positions')
plt.plot(filtered_x[:,1], 'g', label='Estimated Velocity')
plt.title('Positions of tolloco')
plt.xlabel('time step')
plt.ylabel('position')
plt.legend(loc='best')

描画結果 f:id:ray2480:20181118205016p:plain

なんか思ったほどうまくできてない。説明などが色々必要なのだが、後ほど更新する。

2018-11-18

Piaza 始めてみた　〜スキルチェック問題〜

python 競技プログラミング

以前競技プログラミングを少しやってみたいと言ったが、昔Paizaに登録していたのを思い出したのでやってみた。その時のメモ書き

問題

D109:ゾロ目の日付

月と日にちを入力して、ゾロ目だったらYes, そうでなければNoを出力するプログラムを作る。

自分の解答：

md = input().rstrip().split()

m = int(md[0])
d = int(md[1])
m10 = m // 10 
d10 = d // 10
m1 = m % 10 
d1 = d % 10

if(m10 == 0 and d10 == 0 and m1 == d1):
    print("Yes")
elif(m10 == 0 and m1 == d10 and d10 == d1):
    print("Yes")    
elif(d10 == 0 and m10 == m1 and m1 == d1):
    print("Yes")    
elif(m10 == m1 and m1 == d10 and d10 == d1):
    print("Yes") 
else:
    print("No")

うーんこの無理矢理感。２数の１の位、１０の位を取ってきて全部一致するか確認する。２桁目が０（存在しない）ものも考慮するってやり方だけど、桁数増えたらこれどうしようもないでしょみたいな解答。

正規表現使ってやるのが良いんだろうきっと。

pythonでキリ番かどうか判定する - BlankTar

また書き直して作ってみようと思う。

2018-11-18

Tex使い方メモ１

TeX 基礎知識

アクセント記号

アクセント記号: LaTeX

文字頭に打つ記号。ドットとかハットとか。

記号一覧

LaTeXƒRƒ}ƒ“ƒhƒV[ƒgˆê—— - ‚·‚×‚Ä‚Ì‹L†

基本的な記号はだいたい載ってるはず

行列

数式: 行列、ベクトル matrix, vector TeX

\begin{bmatrix}など

LaTeXコマンド - 行列と行列式 - matri, pmatrix, bmatrix

split

数式: split 環境: LaTeX

式を複数行に分ける。

2018-11-18

python基本知識のメモ1

python 基礎知識

python基本知識のメモ1

文字列から数値変換、型確認

print(int('100'))
print(type(int('100')))
# 100
# <class 'int'>

floatでも同じ
Pythonで数字の文字列strを数値int, floatに変換 | note.nkmk.me

型確認はtype(object)で確認。リストや配列なども教えてくれる。
【Python入門】型をtype関数で確認してみよう | 侍エンジニア塾ブログ | プログラミング入門者向け学習情報サイト

ぜろといち

ポンコツ理系大学院生による雑多なブログ

毎日やりたいタスクメモ

英語での住所の書き方

カランメソッドについて情報収拾した事まとめ

カランメソッドの復習方法

カランメソッドは土台作り

カランメソッドの効果とか注意とか

カランレッスン体験談

カルマンフィルタ実装してみた

Piaza 始めてみた　〜スキルチェック問題〜

問題

Tex使い方メモ１

アクセント記号

記号一覧

行列

split

python基本知識のメモ1

if文の「かつ」「または」

切捨て除算

絶対値