20191024のPythonに関する記事は30件です。

投稿日：2019-10-24T23:20:01+09:00

Jupyter Notebook で "NameError: name 'file' is not defined" となる

疑問

Jupyter Notebookで以下のコードを実行してみると...
print(os.path.dirname(__file__))
エラーがでる。
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-8-755a91d3a9c3> in <module>()
      3 > NameError: name '__file__' is not defined
      4 '''
----> 5 print(os.path.dirname(__file__))

NameError: name '__file__' is not defined
どうしてだ...？

原因

JupyterではIPython(REPL)からコードを実行しているが、この環境下だと特殊グローバル変数 __file__は設定されないらしい。

解決

代わりに以下のコードで現在のipynbの絶対パスを取得する。
os.path.abspath('')
参考

stack overflowの投稿

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T21:47:18+09:00

DeepChemでCross-Validationをやってみる

はじめに

DeepChemのサンプルでは、Hold-Outの事例がほとんどで、Cross-Validationの事例が少ない。(DeepChemを用いたベンチマークであるMoleculeNetもHold-Outによる評価となっている）。しかし、Cross-Validationの方が、検証に多くのデータが使えることもあり、高々1000件程度のデータであれば、迷わずCross-Valdationを使いたい。今回、DeepChemでCross-Validationをする方法を調べたのでメモっておく

環境

python 3.6

deepchem 2.2.1.dev54

rdkit 2019.03.3.0

ソースコード

以下がソースコードだ。Splitterクラスのk_fold_splitメソッドにより分割が得られるので、後はscikit-learnのKFoldと同じようにすればよい。
DeelChemCrossValExample.py
import argparse
import numpy as np
import deepchem as dc
from sklearn.ensemble import RandomForestRegressor

def get_model():
    estimator = RandomForestRegressor()
    return dc.models.SklearnModel(estimator, None)

def main():

    parser = argparse.ArgumentParser()
    parser.add_argument("-train", type=str, required=True, help="trainig data file(csv)")
    parser.add_argument("-target_col", type=str, required=True)
    parser.add_argument("-smiles_col", type=str, default="smiles")
    parser.add_argument("-id_col", type=str, default="id")
    parser.add_argument("-cv", type=int, default=5)

    args = parser.parse_args()

    featurizer = dc.feat.RDKitDescriptors()

    # 学習データの読み込み
    loader = dc.data.CSVLoader(tasks=[args.target_col],
                               id_field=args.id_col,
                               smiles_field=args.smiles_col,
                               featurizer=featurizer)

    dataset = loader.featurize(args.train)
    splitter = dc.splits.RandomSplitter(dataset)

    transformers = [dc.trans.NormalizationTransformer(transform_y=True, dataset=dataset)]

    split_datas = splitter.k_fold_split(dataset, args.cv)
    metric = dc.metrics.Metric(dc.metrics.r2_score, np.mean)
    train_scores = []
    validation_scores = []

    # Cross-Validationの実行
    for train_set, validation_set in split_datas:
        model = get_model()
        model.fit(train_set)
        train_score = model.evaluate(train_set, [metric], transformers)
        train_scores.append(train_score)
        validation_score = model.evaluate(validation_set, [metric], transformers)
        validation_scores.append(validation_score)
        print("score train:{0}, val:{1}".format(train_score, validation_score))


if __name__ == "__main__":
    main()
おわりに

手元のscikit-learn版のプログラムと同じデータ、同じアルゴリズムで比較したところほぼ同程度の精度となった。
今回はscikit-learnのRandomForestRegressorをラップしたモデルでやってみたが、DeepChemのアルゴリズムでも同様にできるはずだ。
今後はCross-Validationによるハイパーパラメータ探索なども試してみる予定である。
さて今回でDeepChemを使う準備は概ね整った。次回からは、本格的にDeepChemならではのアルゴリズムを使っていきたい。

参考

-How to implement cross-validation for training set in deepchem? #736

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T21:39:10+09:00

Zendeskのチケット検索で楽するためのツール simple-zendesk-searcher（仮称）を作った

結論

Zendeskのチケットを検索するための簡易的なCLI "simple-zendesk-searcher(仮称)" を作りました。ソースは以下。

https://github.com/hassaku63/simple_zendesk_searcher

今のところはGitHubオンリーでの公開です。現時点でちゃんとしたテストも書けてないし、機能的にもまだまだpypiに登録するほどのツールではないかなと思ったので。

日付検索、タグ検索、リクエスタ検索あたりの機能と、あといい感じのパッケージ名が欲しい。。あと ./output.csv に出力する仕様、これデフォルトはstdoutが良さそうなので近々改変しようと思います。

なぜ作ったのか

自社のビジネスでサービスデスク的なサービスを提供しており、サポートデスクとしてZendeskを採用しています。

AWSのメンテナンス系通知とか、多数のアカウントに一様にばらまかれるタイプの通知が来たとき、とりあえずチケット管理システムには起票されるような実装をすると思います。

で、サポートサービスとしては、その通知に対する「当社見解」的なものを添えて顧客にアナウンスしたかったりします。通知のやり方としては、生成済みチケットに対するリプライの形式がベストでしょう。となると、Zendesk的には複数のチケットに対して一斉に「コメントを送付」するような処理が必要になります。

一斉起票の仕組み自体は自社で既に持っているのですが、対象チケットの抽出部分はZendeskのWeb-UIしか方法がなく、大量のチケットリストを抽出するには些か非効率でした。

で、それならCLIツールとして簡単にでも検索できる仕組みを作ればいいじゃないかと思い、simple-zendesk-searcherを製作しました。

※性質上、ZendeskのAPI-KEYを所持していることがツール利用の必要条件ですが、実際にはAPI-KEY持ってないし持ちたくもない、という立場の方もいるので業務にどう組み込むかの課題は残ってます。そのへんはまた別途考えたい。

使い方 & ソース

GitHubに公開してるので、そちらをご覧ください。

https://github.com/hassaku63/simple_zendesk_searcher

今はキーワード検索のみ対応してます。また、Zendeskの仕様上、スペースを含む文字列で検索したい場合はダブルクォートで囲った文字列を投げる必要があります。お気をつけください。シェルから実行する場合はシングルクォートで外を囲うなりエスケープするなりして、然るべき措置が必要です。

出力するカラムも、なんかこう、カスタマイズの余地を残せる仕様にしたいなぁとは思います。もしくは可能な限りのカラムを出力するような仕様にするか。

まとめ

とりあえずの最低要件を満たすツールは作れたので、あとは現場で使ってもらえるようにAPI-KEY所持してる人に共有していこうかと。とりあえず、日付の条件指定くらいはさくっと実装しよう...

あとこれ、記事の主題とは無関係なんですが。業務遂行の直接的な手段としてのツールを作ることと、それを業務フローに組み込むこと、この2つの間にはかなりの隔たりがあるんですよね。なので、個々のツールも大事とは思いつつも「業務フローに組み込む」ことに対する汎用的な回答をいい感じに作れないかなぁと感じました。

BPM(Business Process Management)システムとの統合を意識するならAPI Gatewayにロジック部分を載せて、実行権限の認可はBPM側に任せるとかが良さそうです。あとは、ツールの実装を変えずに無理なくAPI Gatewayにインテグレーションできるような設計とかできたら幸せになれそうだなと思います。

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T20:52:21+09:00

[ゼロから作るDeep Learning]勾配式実装準備 Pythonで微分と偏微分を実装してみた。

はじめに

この記事はゼロから作るディープラーニング 5章ニューラルネットワークの学習を自分なりに理解して分かりやすくアウトプットしたものです。
文系の自分でも理解することが出来たので、気持ちを楽にして読んでいただけたら幸いです。
また、本書を学習する際に参考にしていただけたらもっと嬉しいです。

微分の知識が必要な理由

ニューラルネットワークの学習では、ニューラルネットワークでの処理から損失関数までを一つの関数と捉えて、各重みやバイアスなどのパラメータを変数と捉えて、各パラメータで偏微分を行うことで最適なパラメータへの方向を導き、少しずつパラメータを最適化していきます。

つまり、ニューラルネットワークの学習において微分の知識は最重要だと言うことです。
今回は、各パラメータの偏微分をまとめたものである勾配を導き出す式勾配式を実装するために、まずは微分と偏微分をPythonで実装したいと思います。

微分

微分とは簡単に説明すると傾きのことです。
微分は関数の変数の値を極少量増加させた時に関数全体の値がどれだけ変化したかを表します。
中高で習った関数の傾きとは変数Xが１増加した時に関数全体の値Yがどれだけ変化したかを表したものでした。
そう考えると、Xの増加量が１から極少量に変わっただけなのでほとんど同じです。

では実装してみます。
# 良い微分の実装
def slopeing_pro(f, x, h = 1e-4):
    return (f(x + h) - f(x - h)) / (2 * h) # 数値微分の誤差を小さくする
上の実装で注目するべき部分は、1e-4 と (f(x + h) - f(x - h)) / (2 * h)です。

1e-4は10の-4乗(0.0001)を表します。この値を極少量として変数に足します。ここで注意点ですが、この値を極端に小さくしすぎてはいけません。(1e-50など)
コンピュータはものすごく小さな数の処理が苦手です。そのため、極端に小さな数を使ってしまうとうまく処理ができなくなってしまいます。

(f(x + h) - f(x - h)) / (2 * h)は真の微分と数値微分との誤差を小さくするために、本来の微分の数式を改良したものです。
(f(x + h) - f(x)) / h のままだと真の微分(瞬間変化量)と数値微分(変数が極少量増加した時の全体変化量)の誤差が大きくなってしまいます。その問題を解消するために上のような改良を施します。

偏微分

偏微分とは、２つ以上の変数が使われた関数(多変数関数)を各変数一つ一つで微分することを言います。
簡単にいうと、一つの変数に絞り、絞った変数にだけ極少量の増加させて微分を求めます。
ニューラルネットワークの学習では、各パラメータを変数とおいて各パラメータで偏微分をして、それをまとめた勾配というものを利用するので、偏微分はすごく重要です。

では実装します。
# 偏微分実装
def slopeings_one(f, x, y, h = 1e-4):
    return (f(x + h, y) - f(x - h, y)) / (2 * h)

# 二変数関数
def two_func(x,y):
    return (4 * x) + (3 * y)

slopeings_one(two_func, 2, 3)
上の実装は二変数関数のxだけに対応した偏微分です。
実際の勾配式には、全ての多変数関数かつその全ての変数に対応できるような実装が求められます。

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T19:55:12+09:00

事前分布に依存しないベイズ推定＆モデル比較

概要・導入

ベイズ統計と頻度主義に基づく統計とを比較したとき、ベイズ統計の欠点として、事後分布が事前分布（prior）に依存することが挙げられる。
これは特に標本数が少ないとき（尤度関数のピークが弱いとき）に顕著になる。
この記事では、arXiv:1910.06646の論文で紹介されている、ベイズ統計における事前分布に依存しないモデル比較の方法についてまとめる。

理論

検証したいモデルを $\mathcal{M}$ ,モデルのパラメータの組を $\Theta=(\theta_1,\theta_2,\dots)$ とおく。
特に、$\Theta$を興味があるパラメータの組 $x$ とそうでないもの（nuisance パラメータ） $\psi$ に分ける：$\Theta\to (x,\psi)$
与えられたデータ $d$ に対し、モデル $\mathcal{M}$ の尤度関数は

\newcommand{\calM}{\mathcal{M}}
\newcommand{\calL}{\mathcal{L}}
{\calL}_{\calM}(\Theta) = {\calL}_{\calM}(x,\psi) = p(d|x,\psi,\calM)

と与えられているとする。

今、 $x$ に対して何かデータ $d$ を用いた信頼区間 or 制限を与えることを考える。
特に、$x=x_0$ なのかそうでないかを考えたいような場合、単純に考えると、これは $x=x_0$ と制限したモデル（ $\mathcal{M}_{x_0}$ とする）と $\mathcal{M}$ とのベイズファクター $B$ （Evidence の比）を見れば良い：

\begin{align*}
B = \frac{Z}{Z_{x_0}} &= \frac{\int d\Theta\ {\calL}_{\calM}(\Theta)\ \pi(\Theta|\calM)}{\int d\psi\ \calL_{\calM_{x_0}}(\psi)\ \pi(\psi|\calM_{x_0})}
\end{align*}

ただしここで

\calL_{\calM_{x_0}}(\psi) \equiv \mathcal{L}_{\calM}(x_0,\psi)

とした。

しかし、これは知りたいパラメータである $x$ についての事前分布のとり方（ $\pi(\Theta|\calM)$ に含まれている）に依存する。
そこで、 $Z$ についても特に $x$ をある値に制限することを考え（このモデルを$\mathcal{M}_x$ とする）、この場合のベイズファクターを考えることにする：

\mathcal{R}(x,x_0|d) \equiv B_{x,x_0} = \frac{Z_x}{Z_{x_0}} = \frac{\int d\psi\ \mathcal{L}(x,\psi)\ \pi(\psi|\calM_x)}{\int d\psi\ \mathcal{L}(x_0,\psi)\ \pi(\psi|\calM_{x_0})}

この時 $\mathcal{R}(x,x_0|d)$ を "relative belief updating ratio"（相対信頼度更新比？）や "shape distortion function"（形状歪曲関数？）と呼ぶ。
$\mathcal{R}(x,x_0|d)$ は $x$ に関する事前分布なしに定まっているので、この意味で事前分布依存性がない（もちろん他のパラメータに対しての事前分布依存性はある）。

実際に$\mathcal{R}(x,x_0|d)$ の値を評価する場合は、Evidence $Z_x,Z_{x_0}$ を計算する代わりに、
ベイズファクターの満たす関係式¹

\frac{p(\mathcal{M_x}|d)}{p({{\calM}_{x_0}}|d)} = B_{x,x_0}\frac{\pi(\mathcal{M_{x}})}{\pi(\calM_{x_0})}

及び

\begin{align}
p({\calM}_x|d) &= p(x|d,\calM)\ p(\calM|d)\\
\pi({\calM}_x) &= \pi(x|\calM)\ \pi(\calM)
\end{align}

と $x$ によらない部分をくくり出せることを用いて導ける

\mathcal{R}(x,x_0|d) = \frac{\frac{p(x|d,M)}{\pi(x|\calM)}}{\frac{p(x_0|d,\calM)}{\pi(x_0|\calM)}}

という表式が便利。

例１：ガウス分布の混合比推定

例として、でかいガウス分布のバックグラウンドにちょこっとだけシグナル成分のガウシアンが乗ったようなデータを考える：

p(x|s) \sim s \mathcal{G}[x;\mu=2,\sigma=0.1] + (1-s) \mathcal{G}[x;\mu=0,\sigma=1]

実際にプロットしてみるとこんな感じ：

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.stats import norm
from scipy.special import logsumexp

n_tot = 1000
n_sig = 50
n_bg = n_tot - n_sig

s_th = n_sig / n_tot

loc_sig = 2
scale_sig = 0.1

data_bg = norm.rvs(size=n_bg)
data_sig = norm.rvs(loc=loc_sig,scale=scale_sig,size=n_sig)
data = pd.Series(np.concatenate([data_bg,data_sig]))

data.hist(bins=64)
plt.show()

このような分布に対して、異なる事前分布（flat, logflat, halfcauchy）を用いて混合比 s_th を推定してみる。
推定に当たっては、EnsembleSampler（単純なMetropolis-Hasting法と比べてステップ幅をいじらなくてもいい感じにサンプリングができる）が便利なemceeを使ってみる。
s_th=0.05くらいならばどちらの事前分布を使ってもそこまで違いはない：

from emcee import EnsembleSampler
import sys

#epsilon = sys.float_info.epsilon
epsilon = 0


#### MCMC hyperparameters ####
nwalkers = 64
dim = 1
threads = 32
n_burnin_repeat = 3
n_burnin = 500
n_mcmc = 2000


#### configration of the prior ####
#loc_halfcauchy = 0
#scale_halfcauchy = 1


def lnlike(param):
    s = param
    p_bg = norm.pdf(data)
    p_sig = norm.pdf(data,loc=loc_sig,scale=scale_sig)
    p_tot = s * p_sig + (1-s) * p_bg
    return np.sum(np.log(p_tot))


def lnprior_flat(param):
    s = param
    if (s < epsilon) or (s > 1):
        return -np.inf
    else:
        return 0


def lnprior_logflat(param):
    s = param
    if (s < epsilon) or (s > 1):
        return -np.inf
    else:
        return -np.log(s)  # prior(s) = 1/x


def lnprior_halfcauchy(param):
    s = param
    if (s < epsilon) or (s > 1):
        return -np.inf
    else:
        return halfcauchy.logpdf(s)


def lnpost_flat(param):
    lnp = lnprior_flat(param)
    if lnp > -np.inf:
        return lnlike(param) + lnp
    else:
        return lnp


def lnpost_logflat(param):
    lnp = lnprior_logflat(param)
    if lnp > -np.inf:
        return lnlike(param) + lnp
    else:
        return lnp


def lnpost_halfcauchy(param):
    lnp = lnprior_halfcauchy(param)
    if lnp > -np.inf:
        return lnlike(param) + lnp
    else:
        return lnp


sampler_flat = EnsembleSampler(
    nwalkers = nwalkers,
    dim = 1,
    lnpostfn = lnpost_flat,
    threads=threads
)

sampler_logflat = EnsembleSampler(
    nwalkers = nwalkers,
    dim = 1,
    lnpostfn = lnpost_logflat,
    threads = threads
)

sampler_halfcauchy = EnsembleSampler(
    nwalkers = nwalkers,
    dim = 1,
    lnpostfn = lnpost_halfcauchy,
    threads = threads
)


pos0_flat = norm.rvs(loc=s_th,scale=0.001,size=nwalkers)[:,np.newaxis]
for i in range(n_burnin_repeat):
    pos0_flat, _, _ = sampler_flat.run_mcmc(pos0_flat,n_burnin)
    sampler_flat.reset()
sampler_flat.run_mcmc(pos0_flat,n_burnin)
res_flat = pd.Series(sampler_flat.flatchain[:,0])


pos0_logflat = norm.rvs(loc=s_th,scale=0.001,size=nwalkers)[:,np.newaxis]
for i in range(n_burnin_repeat):
    pos0_logflat, _, _ = sampler_logflat.run_mcmc(pos0_logflat,n_burnin)
    sampler_logflat.reset()
sampler_logflat.run_mcmc(pos0_logflat,n_burnin)
res_logflat = pd.Series(sampler_logflat.flatchain[:,0])


pos0_halfcauchy = norm.rvs(loc=s_th,scale=0.001,size=nwalkers)[:,np.newaxis]
for i in range(n_burnin_repeat):
    pos0_halfcauchy, _, _ = sampler_halfcauchy.run_mcmc(pos0_halfcauchy,n_burnin)
    sampler_halfcauchy.reset()
sampler_halfcauchy.run_mcmc(pos0_halfcauchy,n_burnin)
res_halfcauchy = pd.Series(sampler_halfcauchy.flatchain[:,0])

出力：s_th=0.05

ress = np.array([res_flat,res_logflat,res_halfcauchy])
bins = np.logspace(np.log10(1e-2),np.log10(ress.max()),256)
res_flat.hist(bins=bins,histtype="step")
res_logflat.hist(bins=bins,histtype="step")
res_halfcauchy.hist(bins=bins,histtype="step")

plt.legend(["flat","logflat","halfcauchy"])
plt.xscale("log")
plt.yscale("log")
plt.show()

事後分布の推定結果：

このくらいだとどれもそんなに変わらない。

しかし、n_tot = 1000, n_sig = 10 (s_th = 0.01)くらいにすると、事前分布の依存性が出てくる：
出力：s_th=0.01

拡大図

さてこの時、 $\mathcal{R}(x,x_0|d)$ を求めてみると：

hist_flat = np.histogram(res_flat,bins,density=True)  # hist(n), bins_edges(n+1)
hist_logflat = np.histogram(res_logflat,bins,density=True)
hist_halfcauchy = np.histogram(res_halfcauchy,bins,density=True)

prior_flat = lambda x: 1
prior_logflat = lambda x: 1/x
prior_halfcauchy = lambda x: halfcauchy.pdf(x,loc=0,scale=scale_halfcauchy)

# relative belief updating ratio (RBUR)
def rbur(hist,prior_func):
    posterior,bins = hist
    r = 0.5
    x = np.exp(r*np.log(bins[1:]) + (1-r)*np.log(bins[:-1]) )
    post_prior = posterior/prior_func(x)
    return post_prior/post_prior[0],x

rbur_flat,x_flat = rbur(hist_flat,prior_flat)
rbur_logflat,_x = rbur(hist_logflat,prior_logflat)
rbur_halfcauchy,_x = rbur(hist_halfcauchy,prior_halfcauchy)

plt.plot(x_flat,rbur_flat)
plt.plot(_x,rbur_logflat)
plt.plot(_x,rbur_halfcauchy)
plt.legend(["flat","logflat","halfcauchy"])
for i in range(1,8,2):
    plt.plot(_x,np.ones_like(_x)*np.exp(-i),"--",c="gray",linewidth=1)
plt.xscale("log")
plt.yscale("log")

plt.xlabel("s")
plt.ylabel(r"$\mathcal{R}$")

こうなる：

確かに $\mathcal{R}$ には事前分布の依存性がなさそう。
もっと頑張ってs_th=0.001にしてみるとこんな感じ：

事前分布の違いによってs==0付近での統計が足らずに若干ふらつくが、それでも十分合っている。

例２：ガウス分布の混合比＆分散推定

上の分布で、更に分布の分散も局外パラメータ (nuisance parameter) として振ってみることにする。
分散の事前分布はいずれも半コーシー分布で固定しとく。
今回もn_sig = 1 としておく。

サンプリング用コード

一次元の時とほとんど同じ。

from emcee import EnsembleSampler
import sys

#epsilon = sys.float_info.epsilon
epsilon = 0


#### MCMC hyperparameters ####
nwalkers = 64
dim = 1
threads = 32
n_burnin_repeat = 3
n_burnin = 500
n_mcmc = 2000


#### configration of the prior ####
#loc_halfcauchy = 0
scale_halfcauchy = 0.5


def lnlike_2d(param):
    s,scale = param
    p_bg = norm.pdf(data)
    p_sig = norm.pdf(data,loc=loc_sig,scale=scale)
    p_tot = s * p_sig + (1-s) * p_bg
    return np.sum(np.log(p_tot))


def lnpost_flat_2d(param):
    s,scale = param
    lnp = lnprior_flat(s) + lnprior_halfcauchy(scale) 
    if lnp > -np.inf:
        return lnlike_2d(param) + lnp
    else:
        return lnp


def lnpost_logflat_2d(param):
    s,scale = param
    lnp = lnprior_logflat(s) + lnprior_halfcauchy(scale)
    if lnp > -np.inf:
        return lnlike_2d(param) + lnp
    else:
        return lnp


def lnpost_halfcauchy_2d(param):
    s,scale = param
    lnp = lnprior_halfcauchy(s) + lnprior_halfcauchy(scale)
    if lnp > -np.inf:
        return lnlike_2d(param) + lnp
    else:
        return lnp


sampler_flat_2d = EnsembleSampler(
    nwalkers = nwalkers,
    dim = 2,
    lnpostfn = lnpost_flat_2d,
    threads=threads
)

sampler_logflat_2d = EnsembleSampler(
    nwalkers = nwalkers,
    dim = 2,
    lnpostfn = lnpost_logflat_2d,
    threads = threads
)

sampler_halfcauchy_2d = EnsembleSampler(
    nwalkers = nwalkers,
    dim = 2,
    lnpostfn = lnpost_halfcauchy_2d,
    threads = threads
)


pos0_flat_2d = norm.rvs(loc=s_th,scale=0.001,size=(nwalkers,2))
for i in range(n_burnin_repeat):
    pos0_flat_2d, _, _ = sampler_flat_2d.run_mcmc(pos0_flat_2d,n_burnin)
    sampler_flat_2d.reset()
sampler_flat_2d.run_mcmc(pos0_flat_2d,n_burnin)
res_flat_2d = pd.Series(sampler_flat_2d.flatchain[:,0])


pos0_logflat_2d = norm.rvs(loc=s_th,scale=0.001,size=(nwalkers,2))
for i in range(n_burnin_repeat):
    pos0_logflat_2d, _, _ = sampler_logflat_2d.run_mcmc(pos0_logflat_2d,n_burnin)
    sampler_logflat_2d.reset()
sampler_logflat_2d.run_mcmc(pos0_logflat_2d,n_burnin)
res_logflat_2d = pd.Series(sampler_logflat_2d.flatchain[:,0])


pos0_halfcauchy_2d = norm.rvs(loc=s_th,scale=0.001,size=(nwalkers,2))
for i in range(n_burnin_repeat):
    pos0_halfcauchy_2d, _, _ = sampler_halfcauchy_2d.run_mcmc(pos0_halfcauchy_2d,n_burnin)
    sampler_halfcauchy_2d.reset()
sampler_halfcauchy_2d.run_mcmc(pos0_halfcauchy_2d,n_burnin)
res_halfcauchy_2d = pd.Series(sampler_halfcauchy_2d.flatchain[:,0])

結果のプロット：

一次元の場合と同じようにプロットしてみる。ただし、局外パラメータの情報は落として、$p(s|d)$ だけ見る。

ress_2d = np.array([res_flat_2d,res_logflat_2d,res_halfcauchy_2d])
bins = np.logspace(np.log10(1e-4),np.log10(ress_2d.max()),64)
res_flat_2d.hist(bins=bins,histtype="step")
res_logflat_2d.hist(bins=bins,histtype="step")
res_halfcauchy_2d.hist(bins=bins,histtype="step")

plt.legend(["flat","logflat","halfcauchy"],loc="upper left")
plt.xscale("log")
plt.yscale("log")
plt.show()

散布図：

局外パラメータの振る舞いについても気になるので一応見ておく：

for sampler in (sampler_flat_2d,sampler_logflat_2d,sampler_halfcauchy_2d):
    plt.xscale("log")
    plt.yscale("log")
    plt.xlabel("s")
    plt.ylabel("scale")
    plt.scatter(*sampler.flatchain.T,s=1,linewidths=0,c=sampler.flatlnprobability)
    plt.colorbar()
    plt.show()

Ralative belief updating function:

問題のやつ。局外パラメータが合っても確かに事前分布依存性が取り除けているか？

hist_flat_2d = np.histogram(res_flat_2d,bins,density=True)  # hist(n), bins_edges(n+1)
hist_logflat_2d = np.histogram(res_logflat_2d,bins,density=True)
hist_halfcauchy_2d = np.histogram(res_halfcauchy_2d,bins,density=True)

prior_flat = lambda x: 1
prior_logflat = lambda x: 1/x
prior_halfcauchy = lambda x: halfcauchy.pdf(x,loc=0,scale=scale_halfcauchy)

# relative belief updating ratio (RBUR)
def rbur(hist,prior_func):
    posterior,bins = hist
    r = 0.5
    x = np.exp(r*np.log(bins[1:]) + (1-r)*np.log(bins[:-1]) )
    post_prior = posterior/prior_func(x)
    return x,post_prior/post_prior[0]

rbur_flat_2d = rbur(hist_flat_2d,prior_flat)
rbur_logflat_2d = rbur(hist_logflat_2d,prior_logflat)
rbur_halfcauchy_2d = rbur(hist_halfcauchy_2d,prior_halfcauchy)

plt.plot(*rbur_flat_2d)
plt.plot(*rbur_logflat_2d)
plt.plot(*rbur_halfcauchy_2d)
plt.legend(["flat","logflat","halfcauchy"])
for i in range(1,8,2):
    plt.plot(rbur_flat_2d[0],np.ones_like(rbur_flat_2d[0])*np.exp(-i),"--",c="gray",linewidth=1)
plt.xscale("log")
plt.yscale("log")

plt.xlabel("s")
plt.ylabel(r"$\mathcal{R}$")

結果：

やはり低サンプリング数に伴う誤差はあるが、確かに事前分布に依存しない様子が見て取れる。

結局なんなのよ

結局のところ $\mathcal{R}$ を使って頻度主義統計（Frequentist）における尤度比検定のベイジアン的な一般化をしてるだけ。

まとめ

尤度比もどきを使うことで事前分布によらない？ベイジアン的な解析ができる。

ただし、ここで $x,\,\psi$ が互いに独立だと仮定（ $\pi(x,\psi|\calM)=\pi(x|\calM)\ \pi(x|\calM)$ ）していることに注意。互いに独立でない場合はこの式は成り立たない。 ↩

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T18:20:20+09:00

pylint をインストールしようとして sitecustomize のエラーに悩まされた話

概要

Visual Studio Code で Python 拡張機能をインストールしたら

Linter pylint is not installed.

と出たので，指示通り pylint をインストールしようとしたら sitecustomize 絡みのエラーが出た．エラーを解消するには usercustomize.py を自分で書いて設置すればよい．

なお，筆者の Python バージョンは下記の通り．
Python 3.8.0 (tags/v3.8.0:fa919fd, Oct 14 2019, 19:37:50) [MSC v.1916 64 bit (AMD64)] on win32
詳細

Visual Studio Code のダイアログに従って pylint をインストールしようとすると
>python -m pip install -U pylint --user
と同等のコマンドが実行される．ところが，Preparing wheel metadata ... のところで次のようなエラーを吐いた．
ERROR: Command errored out with exit status 1:
     command: 'C:\Users\ユーザー名\AppData\Local\Programs\Python\Python38\python.exe' 'C:\Users\ユーザー名\AppData\Local\Programs\Python\Python38\lib\site-packages\pip\_vendor\pep517\_in_process.py' prepare_metadata_for_build_wheel 'C:\Temp\tmpv4cd53yk'
         cwd: C:\Temp\pip-install-6tndi160\lazy-object-proxy
    Complete output (12 lines):
    running dist_info
    creating C:\Temp\pip-modern-metadata-pfp85tml\lazy_object_proxy.egg-info
    writing C:\Temp\pip-modern-metadata-pfp85tml\lazy_object_proxy.egg-info\PKG-INFO
    writing top-level names to C:\Temp\pip-modern-metadata-pfp85tml\lazy_object_proxy.egg-info\top_level.txt
    writing manifest file 'C:\Temp\pip-modern-metadata-pfp85tml\lazy_object_proxy.egg-info\SOURCES.txt'
    reading manifest file 'C:\Temp\pip-modern-metadata-pfp85tml\lazy_object_proxy.egg-info\SOURCES.txt'
    writing manifest file 'C:\Temp\pip-modern-metadata-pfp85tml\lazy_object_proxy.egg-info\SOURCES.txt'
    creating 'C:\Temp\pip-modern-metadata-pfp85tml\lazy_object_proxy.dist-info'
    Error in sitecustomize; set PYTHONVERBOSE for traceback:
    SyntaxError: (unicode error) 'utf-8' codec can't decode byte 0xe0 in position 0: unexpected end of data (sitecustomize.py, line 7)      
    error: invalid command 'bdist_wheel'
エラーメッセージを見ると sitecustomize.py に問題がありそうなのだが，このファイルが存在するはずの C:\Users\ユーザー名\AppData\Local\Programs\Python\Python38\Lib\site-packages にファイルそのものが存在しなかった．

色々調べた結果

最近のUnicodeDecodeError - Qiita

の通りに usercustomize.py を作成し，先ほどの C:\Users\ユーザー名\AppData\Local\Programs\Python\Python38\Lib\site-packages に設置したら，無事 pylint のインストールができた．

参考

最近のUnicodeDecodeError - Qiita

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T17:47:37+09:00

Python よくわかるように書いた問題集no.1

問題1

円の半径を入力し、円の面積と円周の長さを求めるプログラムをかけ。

*
*
*
*
*
*
*
*
*
*
*
*
*
*
#問題1

import math
pi = math.pi
radius = float(input())

area = pi * radius ** 2
length = 2 * pi * radius

print('円の面積は{0:.2f}です。'.format(area))
print('円の円周の長さは、{0:.2f}です。'.format(length))
出力結果
10
円の面積は314.16です。
円の円周の長さは、62.83です。
ソースコードの説明

import mathで円周率を使うためにmathモジュールをインポートする。そして、input()関数でstr型のデータを取り出し、float()関数でfloat型に直す。そして、面積、円周の長さを計算して、'{0:.2f}' .format()で桁数を指定しており、小数点第3位を四捨五入している。

変数の定義

変数名意味

pi 円周率

radius 半径

area 面積

length 円周の長さ

フローチャート
graph LR
A(start) -->B[radiusを入力]
B --> C[area/面積/を計算]
C --> D[length/円周の長さ/を計算]
D --> E[計算結果を表示]
E --> F(end)

変数名	意味
pi	円周率
radius	半径
area	面積
length	円周の長さ

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T17:39:53+09:00

大規模日本語ビジネスニュースコーパスを学習したXLNet（MeCab+Sentencepiece利用）モデルの紹介

はじめに

以前、日本語のBERT事前学習済モデルとELMo学習モデルの紹介記事を投稿しましたストックマークの森長です。
モデル公開の記事を多くの皆様に読んでいただき、ありがとうございます。

昨今の自然言語処理界？では、事前学習モデルであるBERTの登場を皮切りに、XLNet、RoBERTa、ALBERTと多数のモデルが提案され、SOTAを競いあい、大いに盛り上がっています！
ですが、最先端のモデルは英語や中国語で事前学習されたモデルが多く、日本語で試すにはハードルがかなり高いと感じています。
そこで、今回はBERT、ELMoに続いて、XLNetの日本語事前学習済モデルを公開いたします。

XLNetとは

XLNetとは、自己符号化ベースであるBERTの以下懸念点を解消するために作られた、自己回帰ベースのモデルです。

BERTの[MASK]トークンは、fine-tuningの時に使用しないためpre-trainingとfine-tuningの学習に差異が生じる

予測対象の単語が複数あった場合、BERTは、自己回帰モデルのように予測対象の単語間の依存関係を考慮できない。

【原文】： New York is A city.

【BERTのMLM】：[MASK] [MASK] is A city. (NewとYorkの依存関係を学習できない)

XLNetの論文やアルゴリズム等の紹介はとても丁寧に解説していただいている記事が多数ありますので、それらの記事を読んでいたただけると良いと思います。私も大変お世話になっています。

XLNetは、通常の自己回帰モデルとは異なり、単語の予測順序の入れ替えを行い、順方向・逆方向両方の依存関係を学習できるようにしたことがとても注目すべき点だと思います。

弊社のXLNet事前学習済モデルの特徴

学習元データ

以前公開しましたBERTの事前学習済モデルと同様に学習用データを日本語ビジネスニュース記事としています。そのため、本モデルは、主にビジネスに関するドメインで効果を発揮します。

tokenizer

弊社では、tokenizerにMeCab + NEologd及びSentencepieceを利用して、文章をトークン化しています。
まずMeCabを利用することで形態素として正しく分割し、その上でSentencepieceを行うことで、語彙のカバー範囲を広げています。

文字コード

XLNetの公式githubでは、正規化手法としてNFKDを採用しています。そのため、公式githubのライブライリをそのまま使用すると、日本語の濁点・半濁点がロストしてしまいます。
日本語では、濁点・半濁点の有無で意味が変わる可能性もありますので、弊社ではNFKD及びNFKCで正規化した2つのモデルを作成しています。

例)
[原文]ストックマーク社が新サービスを発表した
→[NFKDで正規化すると]ストックマーク社か新サーヒスを発表した (公式ライブラリを利用)
→[NFKCで正規化すると]ストックマーク社が新サービスを発表した (公式ライブラリの正規化をNFKCに変更)

事前学習済モデルの詳細

事前学習用データトークナイザー語彙数正規化シーケンス長

日本語ビジネスニュース記事(300万記事) MeCab + NEologd 及び Sentencepiece 32,000語 NFKD(公式の正規化) 128

日本語ビジネスニュース記事(300万記事) MeCab + NEologd 及び Sentencepiece 32,000語 NFKC 128

ダウンロード

公式(TensorFlow)とPyTorch版で利用可能なモデルを公開します。

ダウンロードリンク
モデルの内訳は以下の通りです。(NFKD及びNFKCのディレクトリ構成は同様です)

TensorFLow版

model.ckpt.meta、 model.ckpt.index、 model.ckpt.data-00000-of-00001

config.json

spiece.model (Sentencepieceのモデル)

PyTorch版

pytorch_model.bin

config.json

spiece.model (Sentencepieceのモデル)

利用方法

事前学習済モデルの使用にあたっては、TensorFlow版及びPyTorch版で公開されているサンプルスクリプトで動作しますが、以下のようにスクリプト及びライブラリを修正してください。

入力データ(入力文章)をMeCab+NEologdでわかち書きしてください。

NFKC版を利用する場合は、各ライブラリ(TensorFlow版及びPyTorch版)のnomalizeをNFKDからNFKCに変更してください。(ライブラリのコードの位置はバージョンによって変わってしまうので、grepでNFKDを検索して場所を特定したほうが良いかもです)

XLNetの日本語事前学習済モデルを利用した所感

BERTとの比較

弊社の検証用タスク（記事の分類）においてはBERTを1%程度下回る性能となっていました。
また、自己回帰モデルの特性上、BERTと比較して計算時間も増加しており、利便性も低下しています。
とはいえ、単語間の関係性が重要になるタスク（固有表現抽出等）に対して適応した場合は違った結果が期待できるため、今後はそれらのタスクで検証を行っていきます。

NFKDとNFKCの比較

NFKD版とNFKC版の精度を比較しましたが、精度に大きな差はありませんでした。
濁点・半濁点がロストしても以下の理由で精度が変わらないのでは無いかと考えていますので、今後検証していきます。

が→か、で→て等の助詞でのロスト

助詞のため、タスクの精度向上に大きく寄与しない可能性がある。

サービス→サーヒス等の一般名詞や固有名詞でのロスト

濁点・半濁点がロストしても、ほとんどの単語で、単語としては一意のままなので影響がない。濁点・半濁点がロストして違う意味になる単語ならば影響があると考えられる。

まとめ

検証中の内容が多い中、最後までお読みいただき、ありがとうございました。
BERT、ELMoと同様に公開したXLNet事前学習済モデルが、自然言語処理の盛り上がりに少しでも貢献ができると大変嬉しいです。

事前学習用データ	トークナイザー	語彙数	正規化	シーケンス長
日本語ビジネスニュース記事(300万記事)	MeCab + NEologd 及び Sentencepiece	32,000語	NFKD(公式の正規化)	128
日本語ビジネスニュース記事(300万記事)	MeCab + NEologd 及び Sentencepiece	32,000語	NFKC	128

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T16:56:26+09:00

statsmodels:一般化線形モデル

一般化線形モデル(Generalized linear model)

1972年にネルダーとウェダーバーンによって提唱された、一般化線形モデル(GLM)は、一般的な線形回帰のモデルに正規分布以外の分布を応答変数の誤差項にもつことを可能にした回帰モデルである。リンク関数を通して応答変数との関係を線形とし、各々の観測値の変化の大きさが予測値の関数となることを可能にするモデルである。これにより線形回帰のモデルを一般化している。

通常の線形回帰モデルは、目的変数$Y$、説明変数$X$、かく乱項$\epsilon$の関係を

$$Y=\beta X + \epsilon $$
としてとらえたもので、パラメータ$\beta$の推定に最小二乗法を用いるとかく乱項には平均ゼロ、相互に無相関で、等分散であるという制約が付く。一方、GLMではパラメータの推定に最尤法やベイズの手法が用いられ柔軟性をもたせている。

GLMには線形回帰、ポアソン回帰、ロジスティック回帰などが含まれる。さらにその発展形として相関があったり、クラスターのあるデータを分析する一般化推定方程式、一般化混合モデル、滑らかな関数を用いたノンパラメトリック回帰の一般化加法モデルなどがある。

GLMは3つの要素で構成される。

指数分布族の確率分布

線形の予測因子:$\eta = X \beta$

リンク関数$g$ :$E(Y|X) = \mu = g^{−1}(\eta)$

確率変数 $Y$ が指数分布族であるので、確率密度関数 $f(y)$ は、

$$f(y;\theta,\phi)=\exp ( \frac{y\theta -b(\theta )}{\phi }+c(y,\theta) )$$

の形式で表される。ここで

$\theta$: 正準 (canonical) パラメーター

$\phi$: 分散パラメーター

$b(\theta)$, $c(y,\theta)$: スカラー関数

である。正準パラメーター$\theta$が、滑らかな関数 $g(\theta)$ と、別の確率変数 $X$ の実現値 $x$ を用いて、 $g(\theta) = \beta^T x$ と表せるのが、一般化線形モデルである。$b^{'}(\theta)$, $\phi b^"(\theta)$は平均と分散になる。

たとえば、

回帰モデル bθ ϕ c(y,ϕ) リンク関数

線形回帰 $\theta^2 /2$ $\sigma^2$ $−(y²/\sigma^2 + \log 2\pi \sigma^2)/2$ $\theta$

ロジスティック回帰 $-\log(1 − p)$ 1 0 $\theta$

プロビット回帰 $-\log(1 − p)$ 1 0 $\Psi^{−1} (p)$

$Y$が平均$\theta$, 分散$\sigma^2$の正規分布であれば、$b(\theta) = \theta^2 /2$, $\phi = \sigma^2, c(y,\phi) = −(y²/\sigma^2 + \log 2\pi \sigma^2)/2$であればよい。
リンク関数が$g(\theta) = \theta$ (正準リンクとよぶ)であれば、平均 $\theta$ は $\beta^T x$ となり、通常の線形回帰となる。

また、$Y$ が生起確率 $p$ のベルヌーイ分布（1の確率 $p$, 0の確率$1 − p$ ）であれば、$b(\theta) = −\log(1 − p)$ (ただし $p = e^\theta/(1 + e^\theta) )$, $\phi = 1, c = 0$ であればよい。リンク関数を$g(\theta) = \theta$ とすれば、ロジスティック回帰となる。ここで$\theta=\log (p/q)=\log [p/(1-p)]$

また、リンク関数を$g(\theta) = \Psi^{−1} (p)$ ($\Psi$を標準正規分布の累積分布関数)とすれば、プロビット回帰となる。このように一般化線形モデルは柔軟な構造をもつ。

(注)この部分は統計分布ハンドブック、Generalized linear model wikiではパラメータがベクトルの場合とスカラーの場合で分けて議論をしている。そちらの方が分かりやすいかも。

statsmodelsでは指数分布族はつぎのように表現される。

$$ f(y|\theta,\phi,w)=c(y,\phi,w)\exp(y\theta−b(\theta)\phi w)$$

ここで$w$はまだサポートされていない。確率分布とパラメータの関係はつぎの表を参照：

Dist Domain μ=E v(μ)

Binomial B(n,p) 0,1,…,n np μ$-\mu^2$/n

Poisson P(μ)0,1,…,∞ μ μ

Neg. Binom. NB(μ,α) 0,1,…,∞ μ μ+αμ$^2$

Gaussian/Normal N(μ,$\sigma^2$) (−∞,∞) μ 1

Gamma N(μ,ν) (0,∞) μ μ$^2$

Inv. Gauss. IG(μ,σ2) (0,∞) μ μ$^3$

Tweedie p≥1 depends on p μ μ$^p$

Dist θ(μ) bθ ϕ

Binomial B(n,p) log(p/(1−p)) $n\log(1+e^\theta)$ 1

Poisson P(μ) log(μ) e$^\theta$ 1

Neg. Binom. NB(μ,α) log(αμ/(1+αμ)) −1/αlog(1−αe$^\theta$) 1

Gaussian/NormalN(μ,σ2) μ 1/2θ$^2$ $\sigma^2$

Gamma N(μ,ν) −1μ −log(−θ) 1ν

Inv. Gauss. IG (μ,σ2) −1/2μ$^2$ −2$\sqrt{-\theta}$

Tweedie p≥1 μ1−p1−p α−1α(θα−1)α ϕ

変数、パラメータとcodeの関係はつぎの表を参照：

変数コード

Y,y endog

x exog

β params

μ mu

g 分布族に対するリンク引数

ϕ scale

w まだサポートされていない。 (i.e. w=1)

例：教育の評価に関する分析(二項分布)

Jeff Gill(2000)のStar98の教育関連のデータを取得する。データの詳細は
print(sm.datasets.star98.NOTE)
によって得ることができる。

観測データ数：303(カリフォルニア州のカウンティ)

変数の数：13項目＋8相互作用項

変数名と内容：
- NABOVE   - 数学の成績で国内の中央値を上回る生徒の数
- NBELOW   - 数学の成績で国内の中央値を下回る生徒の数
- LOWINC   - 低所得層の生徒の割合
- PERASIAN - アジアの生徒の割合
- PERBLACK - ブラックの生徒の割合
- PERHISP  - ヒスパニックの生徒の割合
- PERMINTE - マイノリティ教師の割合
- AVYRSEXP - 先生の教育に従事した年数の和を先生の数で割った平均値
- AVSALK   - 総給与予算額をフルタイムの先生の数で割った平均値
- PERSPENK - １生徒あたりの支出額
- PTRATIO  - 先生と生徒の比率
- PCTAF    - UC/CSUの準備コースを受けた生徒の割合
- PCTCHRT  - チャータースクールの割合
- PCTYRRND - year-round スクールの割合
相互作用項
- PERMINTE_AVYRSEXP
- PEMINTE_AVSAL
- AVYRSEXP_AVSAL
- PERSPEN_PTRATIO
- PERSPEN_PCTAF
- PTRATIO_PCTAF
- PERMINTE_AVTRSEXP_AVSAL
- PERSPEN_PTRATIO_PCTAF
データを取得して説明変数に定数を加える。
import numpy as np
import statsmodels.api as sm

#Load the data and add a constant to the exogenous (independent) variables:
data = sm.datasets.star98.load(as_pandas=False)
data.exog = sm.add_constant(data.exog, prepend=False)
データの詳細を得る。
print(sm.datasets.star98.NOTE)
Number of Observations - 303 (counties in California).

Number of Variables - 13 and 8 interaction terms.

Definition of variables names::

    NABOVE   - Total number of students above the national median for the
               math section.
    NBELOW   - Total number of students below the national median for the
               math section.
    LOWINC   - Percentage of low income students
    PERASIAN - Percentage of Asian student
    PERBLACK - Percentage of black students
    PERHISP  - Percentage of Hispanic students
    PERMINTE - Percentage of minority teachers
    AVYRSEXP - Sum of teachers' years in educational service divided by the
            number of teachers.
    AVSALK   - Total salary budget including benefits divided by the number
               of full-time teachers (in thousands)
    PERSPENK - Per-pupil spending (in thousands)
    PTRATIO  - Pupil-teacher ratio.
    PCTAF    - Percentage of students taking UC/CSU prep courses
    PCTCHRT  - Percentage of charter schools
    PCTYRRND - Percentage of year-round schools

    The below variables are interaction terms of the variables defined
    above.

    PERMINTE_AVYRSEXP
    PEMINTE_AVSAL
    AVYRSEXP_AVSAL
    PERSPEN_PTRATIO
    PERSPEN_PCTAF
    PTRATIO_PCTAF
    PERMINTE_AVTRSEXP_AVSAL
    PERSPEN_PTRATIO_PCTAF
#2値の従属変数(Success: NABOVE, Failure: NBELOW):
print(data.endog[:5,:])
[[452. 355.]
 [144.  40.]
 [337. 234.]
 [395. 178.]
 [  8.  57.]]
#すべての独立変数・説明変数
print(data.exog[:2,:])
[[3.43973000e+01 2.32993000e+01 1.42352800e+01 1.14111200e+01
  1.59183700e+01 1.47064600e+01 5.91573200e+01 4.44520700e+00
  2.17102500e+01 5.70327600e+01 0.00000000e+00 2.22222200e+01
  2.34102872e+02 9.41688110e+02 8.69994800e+02 9.65065600e+01
  2.53522420e+02 1.23819550e+03 1.38488985e+04 5.50403520e+03
  1.00000000e+00]
 [1.73650700e+01 2.93283800e+01 8.23489700e+00 9.31488400e+00
  1.36363600e+01 1.60832400e+01 5.95039700e+01 5.26759800e+00
  2.04427800e+01 6.46226400e+01 0.00000000e+00 0.00000000e+00
  2.19316851e+02 8.11417560e+02 9.57016600e+02 1.07684350e+02
  3.40406090e+02 1.32106640e+03 1.30502233e+04 6.95884680e+03
  1.00000000e+00]]
# FIT AND SUMMARY
glm_binom = sm.GLM(data.endog, data.exog, family=sm.families.Binomial())
res = glm_binom.fit()
print(res.summary())
 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:           ['y1', 'y2']   No. Observations:                  303
Model:                            GLM   Df Residuals:                      282
Model Family:                Binomial   Df Model:                           20
Link Function:                  logit   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -2998.6
Date:                Thu, 24 Oct 2019   Deviance:                       4078.8
Time:                        14:30:48   Pearson chi2:                 4.05e+03
No. Iterations:                     5                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
x1            -0.0168      0.000    -38.749      0.000      -0.018      -0.016
x2             0.0099      0.001     16.505      0.000       0.009       0.011
x3            -0.0187      0.001    -25.182      0.000      -0.020      -0.017
x4            -0.0142      0.000    -32.818      0.000      -0.015      -0.013
x5             0.2545      0.030      8.498      0.000       0.196       0.313
x6             0.2407      0.057      4.212      0.000       0.129       0.353
x7             0.0804      0.014      5.775      0.000       0.053       0.108
x8            -1.9522      0.317     -6.162      0.000      -2.573      -1.331
x9            -0.3341      0.061     -5.453      0.000      -0.454      -0.214
x10           -0.1690      0.033     -5.169      0.000      -0.233      -0.105
x11            0.0049      0.001      3.921      0.000       0.002       0.007
x12           -0.0036      0.000    -15.878      0.000      -0.004      -0.003
x13           -0.0141      0.002     -7.391      0.000      -0.018      -0.010
x14           -0.0040      0.000     -8.450      0.000      -0.005      -0.003
x15           -0.0039      0.001     -4.059      0.000      -0.006      -0.002
x16            0.0917      0.015      6.321      0.000       0.063       0.120
x17            0.0490      0.007      6.574      0.000       0.034       0.064
x18            0.0080      0.001      5.362      0.000       0.005       0.011
x19            0.0002   2.99e-05      7.428      0.000       0.000       0.000
x20           -0.0022      0.000     -6.445      0.000      -0.003      -0.002
const          2.9589      1.547      1.913      0.056      -0.073       5.990
==============================================================================
リンク関数はロジット関数、最適化の方法は反復再重み付け最小二乗法 (IRLS法; iteratively reweighted least squares method)が用いられる。これは再重み付けしながら，重み付2乗和誤差を最小化する。Deviance(逸脱度)は予測誤差の平方和のようなもので、それぞれの分布関数により異なる。

$D = 2\sum_i $(freq_weights_i$ * $var_weights$ *(llf(endog_i, endog_i) - llf(endog_i, \mu_i)))$
#分析対象の結果：試行の数(NABOVEの和),パラメータ、t－値
print('Total number of trials:',  data.endog[0].sum())
print('Parameters: ', res.params)
print('T-values: ', res.tvalues)
Total number of trials: 807.0
Parameters:  [-1.68150366e-02  9.92547661e-03 -1.87242148e-02 -1.42385609e-02
  2.54487173e-01  2.40693664e-01  8.04086739e-02 -1.95216050e+00
 -3.34086475e-01 -1.69022168e-01  4.91670212e-03 -3.57996435e-03
 -1.40765648e-02 -4.00499176e-03 -3.90639579e-03  9.17143006e-02
  4.89898381e-02  8.04073890e-03  2.22009503e-04 -2.24924861e-03
  2.95887793e+00]
T-values:  [-38.74908321  16.50473627 -25.1821894  -32.81791308   8.49827113
   4.21247925   5.7749976   -6.16191078  -5.45321673  -5.16865445
   3.92119964 -15.87825999  -7.39093058  -8.44963886  -4.05916246
   6.3210987    6.57434662   5.36229044   7.42806363  -6.44513698
   1.91301155]
#First differences: すべての説明変数の平均値を一定として、
# 反応変数への影響を得るために低所得層の割合を変化させる。

from scipy import stats
means = data.exog.mean(axis=0)
means25 = means.copy()
means25[0] = stats.scoreatpercentile(data.exog[:,0], 25)
means75 = means.copy()
means75[0] = lowinc_75per = stats.scoreatpercentile(data.exog[:,0], 75)
resp_25 = res.predict(means25)
resp_75 = res.predict(means75)
diff = resp_75 - resp_25

#The interquartile first difference for the percentage of 
#low income households in a school district is:
print("%2.4f%%" % (diff*100))

-11.8753%
%matplotlib inline
from matplotlib import pyplot as plt

# Plots p=NABOVE/(NAOVE+NBELO)とその期待値をプロット

nobs = res.nobs
y = data.endog[:,0]/data.endog.sum(1)
yhat = res.mu
from statsmodels.graphics.api import abline_plot
fig, ax = plt.subplots()
ax.scatter(yhat, y)
line_fit = sm.OLS(y, sm.add_constant(yhat, prepend=True)).fit()
abline_plot(model_results=line_fit, ax=ax)


ax.set_title('Model Fit Plot')
ax.set_ylabel('Observed values')
ax.set_xlabel('Fitted values');
fig, ax = plt.subplots()

ax.scatter(yhat, res.resid_pearson)
ax.hlines(0, 0, 1)
ax.set_xlim(0, 1)
ax.set_title('Residual Dependence Plot')
ax.set_ylabel('Pearson Residuals')
ax.set_xlabel('Fitted values')
# Pearson residualsは(endog - mu)/sqrt(VAR(mu))でVARは分布特有の分散
from scipy import stats

fig, ax = plt.subplots()

resid = res.resid_deviance.copy()
resid_std = stats.zscore(resid)
ax.hist(resid_std, bins=25)
ax.set_title('Histogram of standardized deviance residuals');
from statsmodels import graphics
graphics.gofplots.qqplot(resid, line='r')
例：スコットランの選挙 (ガンマ分布)

データはGillの例を基にしている。地方税に関する権限を議会に与えるために「はい」と答えた人の割合である。データは32の地区に分割されている。この例の説明変数には1997年4月に徴収された2人の成人当たりの調整前の地方税の額、1998年1月の失業保険総給付額の女性の割合、標準死亡率(UK=100)、労働組合参加率、地域GDP、5から15歳までの子供の割合、女性の失業と地方税の関係が示されている。Gamma関数は比例計数応答変数に対応している(Gamma for proportional count response)。

観測データ数：32

変数の数：8項目
- YES    - 「はい」と答えた人の割合
- COUTAX - 地方税の額
- UNEMPF - 失業保険総給付額の女性の割合
- MOR    - 標準化死亡率(UK=100)
- ACT    - 労働組合参加率
- GDP    - 地域GDP
- AGE    - 5から15歳までの子供の割合
- COUTAX_FEMALEUNEMP - 女性の失業と地方税の関係
説明変数には定数を付加した。

データの詳細はつぎのように得られる。
print(sm.datasets.scotland.DESCRLONG)
This data is based on the example in Gill and describes the proportion of
voters who voted Yes to grant the Scottish Parliament taxation powers.
The data are divided into 32 council districts. This example's explanatory
variables include the amount of council tax collected in pounds sterling as
of April 1997 per two adults before adjustments, the female percentage of
total claims for unemployment benefits as of January, 1998, the standardized
mortality rate (UK is 100), the percentage of labor force participation,
regional GDP, the percentage of children aged 5 to 15, and an interaction term
between female unemployment and the council tax.

The original source files and variable information are included in
/scotland/src/
data2 = sm.datasets.scotland.load()
data2.exog = sm.add_constant(data2.exog, prepend=False)
print(data2.exog[:5,:])
print(data2.endog[:5])
[[7.12000e+02 2.10000e+01 1.05000e+02 8.24000e+01 1.35660e+04 1.23000e+01
  1.49520e+04 1.00000e+00]
 [6.43000e+02 2.65000e+01 9.70000e+01 8.02000e+01 1.35660e+04 1.53000e+01
  1.70395e+04 1.00000e+00]
 [6.79000e+02 2.83000e+01 1.13000e+02 8.63000e+01 9.61100e+03 1.39000e+01
  1.92157e+04 1.00000e+00]
 [8.01000e+02 2.71000e+01 1.09000e+02 8.04000e+01 9.48300e+03 1.36000e+01
  2.17071e+04 1.00000e+00]
 [7.53000e+02 2.20000e+01 1.15000e+02 6.47000e+01 9.26500e+03 1.46000e+01
  1.65660e+04 1.00000e+00]]
[60.3 52.3 53.4 57.  68.7]
#Fit and summary

data = sm.datasets.scotland.load(as_pandas=False)
data.exog = sm.add_constant(data.exog)
gamma_model = sm.GLM(data.endog, data.exog, family=sm.families.Gamma())
gamma_results = gamma_model.fit()
print(gamma_results.summary())
 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                      y   No. Observations:                   32
Model:                            GLM   Df Residuals:                       24
Model Family:                   Gamma   Df Model:                            7
Link Function:          inverse_power   Scale:                       0.0035843
Method:                          IRLS   Log-Likelihood:                -83.017
Date:                Thu, 24 Oct 2019   Deviance:                     0.087389
Time:                        12:25:13   Pearson chi2:                   0.0860
No. Iterations:                     6                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.0178      0.011     -1.548      0.122      -0.040       0.005
x1          4.962e-05   1.62e-05      3.060      0.002    1.78e-05    8.14e-05
x2             0.0020      0.001      3.824      0.000       0.001       0.003
x3         -7.181e-05   2.71e-05     -2.648      0.008      -0.000   -1.87e-05
x4             0.0001   4.06e-05      2.757      0.006    3.23e-05       0.000
x5         -1.468e-07   1.24e-07     -1.187      0.235   -3.89e-07    9.56e-08
x6            -0.0005      0.000     -2.159      0.031      -0.001   -4.78e-05
x7         -2.427e-06   7.46e-07     -3.253      0.001   -3.89e-06   -9.65e-07
==============================================================================
例：非正準リンクをもつガウス分布

人工的なデータを作成、 family=sm.families.Gaussian(sm.families.links.log)として指定している。
nobs2 = 100
x = np.arange(nobs2)
np.random.seed(54321)
X = np.column_stack((x,x**2))
X = sm.add_constant(X, prepend=False)
lny = np.exp(-(.03*x + .0001*x**2 - 1.0)) + .001 * np.random.rand(nobs2)

gauss_log = sm.GLM(lny, X, family=sm.families.Gaussian(sm.families.links.log))
gauss_log_results = gauss_log.fit()
print(gauss_log_results.summary())
  Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                      y   No. Observations:                  100
Model:                            GLM   Df Residuals:                       97
Model Family:                Gaussian   Df Model:                            2
Link Function:                    log   Scale:                      1.0531e-07
Method:                          IRLS   Log-Likelihood:                 662.92
Date:                Thu, 24 Oct 2019   Deviance:                   1.0215e-05
Time:                        12:26:19   Pearson chi2:                 1.02e-05
No. Iterations:                     7                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
x1            -0.0300    5.6e-06  -5361.316      0.000      -0.030      -0.030
x2         -9.939e-05   1.05e-07   -951.091      0.000   -9.96e-05   -9.92e-05
const          1.0003   5.39e-05   1.86e+04      0.000       1.000       1.000
==============================================================================
関連・参考サイト

statsmodels: 一般化加法モデル
 Welcome to Statsmodels’s Documentation
一般化線形モデルwiki
一般化線形モデルwiki(英語)
2 4.devianceと尤度比検定

参考

「Python3ではじめるシステムトレード」(パンローリング)

「統計分布ハンドブック」(朝倉書店)

セミナー案内

ゼロからまなぶ統計学（Python統計ソフト statsmodels入門）
誰でも簡単統計ソフト
https://kinyuzaimu.connpass.com/event/149178/

回帰モデル	bθ	ϕ	c(y,ϕ)	リンク関数
線形回帰	$\theta^2 /2$	$\sigma^2$	$−(y²/\sigma^2 + \log 2\pi \sigma^2)/2$	$\theta$
ロジスティック回帰	$-\log(1 − p)$	1	0	$\theta$
プロビット回帰	$-\log(1 − p)$	1	0	$\Psi^{−1} (p)$

Dist	Domain	μ=E	v(μ)
Binomial B(n,p)	0,1,…,n	np	μ$-\mu^2$/n
Poisson P(μ)0,1,…,∞	μ	μ
Neg. Binom. NB(μ,α)	0,1,…,∞	μ	μ+αμ$^2$
Gaussian/Normal N(μ,$\sigma^2$)	(−∞,∞)	μ	1
Gamma N(μ,ν)	(0,∞)	μ	μ$^2$
Inv. Gauss. IG(μ,σ2)	(0,∞)	μ	μ$^3$
Tweedie p≥1	depends on p	μ	μ$^p$

Dist	θ(μ)	bθ	ϕ
Binomial B(n,p)	log(p/(1−p))	$n\log(1+e^\theta)$	1
Poisson P(μ)	log(μ)	e$^\theta$	1
Neg. Binom. NB(μ,α)	log(αμ/(1+αμ))	−1/αlog(1−αe$^\theta$)	1
Gaussian/NormalN(μ,σ2)	μ	1/2θ$^2$	$\sigma^2$
Gamma N(μ,ν)	−1μ	−log(−θ)	1ν
Inv. Gauss. IG	(μ,σ2)	−1/2μ$^2$	−2$\sqrt{-\theta}$
Tweedie p≥1	μ1−p1−p	α−1α(θα−1)α	ϕ

変数	コード
Y,y	endog
x	exog
β	params
μ	mu
g	分布族に対するリンク引数
ϕ	scale
w	まだサポートされていない。 (i.e. w=1)

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T16:01:33+09:00

flask_sqlalchemyで複数のデータベースに接続したい

flask_sqlalchemyを使って複数のデータベースに接続する場合の書き方です。

以下はmysqlを３つ接続する場合です。

ディレクトリ構成

├── startup.py  #起動ファイル
├── application
│   ├── __init__.py  #アプリケーション本体ファイル
│   ├── models.py
│   └── views.py
├── instance
│   └── config.py

コンフィグの書き方

１つ目のデータベースをSQLALCHEMY_DATABASE_URIへ記載しています。
残りの２つのデータベースをSQLALCHEMY_BINDSに辞書型で指定しています。

config.py

DB_USER = 'test_user0'
DB_PASS = 'password0'
DB_HOST = 'xxx.xxx.xxx.101'
DB_NAME = 'db1'
db_uri = "mysql+pymysql://{0}:{1}@{2}/{3}?charset=utf8".format(DB_USER, DB_PASS, DB_HOST, DB_NAME)

DB_USER1 = 'test_user1'
DB_PASS1 = 'password1'
DB_HOST1 = 'xxx.xxx.xxx.201'
DB_NAME1 = 'db234'
db_uri_another1 = "mysql+pymysql://{0}:{1}@{2}/{3}?charset=utf8".format(DB_USER1, DB_PASS1, DB_HOST1, DB_NAME1)

DB_USER2 = 'test_user2'
DB_PASS2 = 'password2'
DB_HOST2 = 'xxx.xxx.abb.202'
DB_NAME2 = 'db456'
db_uri_another2 = "mysql+pymysql://{0}:{1}@{2}/{3}?charset=utf8".format(DB_USER2, DB_PASS2, DB_HOST2, DB_NAME2)

SQLALCHEMY_DATABASE_URI = db_uri
SQLALCHEMY_BINDS = {"adding_db1": db_uri_another1, "adding_db2":db_uri_another2}

モデルの書き方

adding_db1のテーブル(table1)を利用したい場合は
SQLALCHEMY_BINDSで指定したbind_keyをtablenameの前に書きます。

models.py

from application import db

class sample01(db.Model):
   __bind_key__ = 'adding_db1'
   __tablename__ = 'table1'
   id = db.Column(db.Integer, primary_key=True, autoincrement=True)
   name = db.Column(db.String(30))
   memo = db.Column(db.String(150))

※bind_keyを書かなかった場合はSQLALCHEMY_DATABASE_URIに接続するようです。

（参考）

startup.py

from application import app


if __name__ == "__main__":
    app.run(host='0.0.0.0', port=5000)

___init__.py

from flask import Flask
from flask_sqlalchemy import SQLAlchemy

app = Flask(__name__, instance_relative_config=True) #コンフィグはinstanceフォルダを見るように指定
app.config.from_pyfile('config.py', silent=False) #config.pyが見つからない場合はエラー
db = SQLAlchemy(app)

import application.views
import application.models

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T15:43:40+09:00

printやdebugやpdbはなんかいや

参考

僕みたいにpython初心者の人に知って欲しい

loopとか回しているときになんでエラー出んねんーーー

クラス作って一気にメソッド大量に作ったけど動くんかこれーーー

ってなる君、君は今からメモするテストについて知るべきだ

　単純なテスト

def astest():
    assert(11 in range(10)), 'this is assert test'

if __name__ == '__main__':
    astest()

出力

Traceback (most recent call last):
  File "/Users/kiwamizamurai/Desktop/uniuni.py", line 41, in <module>
    astest()
  File "/Users/kiwamizamurai/Desktop/uniuni.py", line 32, in astest
    assert(11 in range(10)), 'this is assert test'
AssertionError: this is assert test

assert関数の引数がTrueだと何も問題が起きない
そしてFalseの場合はAssertionErrorが発生する

変数の型のチェックなど単純なものはこれで十分戦える

最強のユニットテスト

もっと大きく、クラス規模でテストしたいときはunittestを使おう

次のメソッドを使う

メソッド	確認事項
assertEqual(a, b)	a == b
assertNotEqual(a, b)	a != b
assertTrue(x)	bool(x) is True
assertFalse(x)	bool(x) is False
assertIs(a, b)	a is b
assertIsNot(a, b)	a is not b
assertIsNone(x)	x is None
assertIsNotNone(x)	x is not None
assertIn(a, b)	a in b
assertNotIn(a, b)	a not in b
assertIsInstance(a, b)	isinstance(a, b)
assertNotIsInstance(a, b)	not isinstance(a, b)

import random
import unittest

class TestSequenceFunctions(unittest.TestCase):

    def setUpClass():
        print('*** 全体前処理 ***')

    def tearDownClass():
        print('*** 全体後処理 ***')

    # INITIALIZE
    def setUp(self):
        # ここに各テスト関数を実行する前に呼び出す処理を書く。
        self.seq = range(10)
        print('\n*** test-前-処理 ***')


    def tearDown(self):
        # ここに各テスト関数を実行した後に呼び出す処理を書く。
        print('\n*** test-後-処理 ***')


    # testで始まるメソッド名を記述する

    # TEST 1
    def test_shuffle(self):
        random.shuffle(self.seq)
        self.seq.sort()
        self.assertEqual(self.seq, range(10))
        self.assertRaises(TypeError, random.shuffle, (1,2,3))

    # TEST 2
    def test_choice(self):
        element = random.choice(self.seq)
        # self.assertTrue(element in self.seq)
        self.assertTrue(element in range(10,20))

    @unittest.skip("skpyはtestをスキップ")
    def test_skipy(self):
        print('skpyメソッドだよ')



if __name__ == '__main__':
    # main()で実行できる。
    # unittest.main()

    # 個別に実行するテストを取得する。
    # test_shuffle = TestSequenceFunctions('test_shuffle')
    # test_choice = TestSequenceFunctions('test_choice')

    # テストスイートに登録してランナーでまとめて実行できる。
    # TestSuite = unittest.TestSuite()
    # TestSuite.addTest(test_shuffle)
    # TestSuite.addTest(test_choice)
    # runner = unittest.TextTestRunner()
    # runner.run(TestSuite)

    # 出力するメッセージレベルはverbosity引数で設定
    suite = unittest.TestLoader().loadTestsFromTestCase(TestSequenceFunctions)
    unittest.TextTestRunner(verbosity=2).run(suite)

TestLoaderの場合の出力

*** 全体前処理 ***
test_choice (__main__.TestSequenceFunctions) ... 
*** test-前-処理 ***

*** test-後-処理 ***
FAIL
test_shuffle (__main__.TestSequenceFunctions) ... 
*** test-前-処理 ***

*** test-後-処理 ***
ERROR
test_skipy (__main__.TestSequenceFunctions) ... skipped 'skpyはtestをスキップ'

======================================================================
ERROR: test_shuffle (__main__.TestSequenceFunctions)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Users/kiwamizamurai/Desktop/uniuni.py", line 28, in test_shuffle
    random.shuffle(self.seq)
  File "/Users/kiwamizamurai/.pyenv/versions/3.7.4/lib/python3.7/random.py", line 278, in shuffle
    x[i], x[j] = x[j], x[i]
TypeError: 'range' object does not support item assignment

======================================================================
FAIL: test_choice (__main__.TestSequenceFunctions)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Users/kiwamizamurai/Desktop/uniuni.py", line 37, in test_choice
    self.assertTrue(element in range(10,20))
AssertionError: False is not true

*** 全体後処理 ***
----------------------------------------------------------------------
Ran 3 tests in 0.001s

FAILED (failures=1, errors=1, skipped=1)

unittest.main()の場合の出力

*** 全体前処理 ***

*** test-前-処理 ***

*** test-後-処理 ***
F
*** test-前-処理 ***

*** test-後-処理 ***
Es
======================================================================
ERROR: test_shuffle (__main__.TestSequenceFunctions)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Users/kiwamizamurai/Desktop/uniuni.py", line 28, in test_shuffle
    random.shuffle(self.seq)
  File "/Users/kiwamizamurai/.pyenv/versions/3.7.4/lib/python3.7/random.py", line 278, in shuffle
    x[i], x[j] = x[j], x[i]
TypeError: 'range' object does not support item assignment

======================================================================
FAIL: test_choice (__main__.TestSequenceFunctions)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Users/kiwamizamurai/Desktop/uniuni.py", line 37, in test_choice
    self.assertTrue(element in range(10,20))
AssertionError: False is not true

----------------------------------------------------------------------
*** 全体後処理 ***
Ran 3 tests in 0.001s

FAILED (failures=1, errors=1, skipped=1)

特定のメソッドのみ実行の場合の出力

self.assertTrue(element in range(1,20))

エラーが出ないように上のように書き換えてから実行すると

*** 全体前処理 ***

*** test-前-処理 ***

*** test-後-処理 ***
*** 全体後処理 ***
.
----------------------------------------------------------------------
Ran 1 test in 0.000s

OK

スキップできるのは良いねええ

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T15:38:59+09:00

保有数が0の場合の売り注文を回避する #QuantX

order is ignored due to amount is 0

QuantXがバージョンアップしてからこんな警告が出るようになりました。

警告文の一部を抜粋
2018-07-12 jp.stock.9984 order is ignored due to amount is 0.
保有数が0で売り注文を出してるから注文を無視するという警告文が出ています。
バックのシステムが無視してくれていますが、黄色い警告文が出ているので不安になります。
handle_signals関数下で処理して、警告文が出ないようにしたいと思います。

警告文が出る条件

色々試行錯誤してみた結果、以下の条件下で出るそうです。

1. 保有数が0で売り注文を出している。

QuantXは現時点ではロングでの入りしかできません（ショートから入ることはできません。）
つまり、保有数が0の時売却する株がないため売ることができないため。システムが自動的に取引を無視してくれます。

2. 株を購入する際、売買条件にあっていない。

例えば、株価2000円の株を単元株数買うと20万円必要になります。
ここで、金額指定を10万円に設定していたりするとシステムが自動的に取引を無視してくれます。

警告文はなるべく消したい。

システムが自動処理してくれるとはいえ、警告文が出るのはやはり嫌なものです。
上記の２に関してはアルゴリズム作成の際の銘柄選択にあった注文方法を選ぶのが早そうですが、1に関してはプログラムで解決できそうなので実装していきます。

set型を使う

今回はPythonのset型を使い、ポジションを持っていない銘柄を判別していきます。
そもそもset型とは集合を扱うための型であり
・重複した要素がない
・要素に順番がない
という二つの特徴を持つ複数の値を格納できる型です。

方針

方法はシンプルです。

1.portfolio_positionsにない銘柄をnone_symsと定義したset型に入れる。
2.portfolio_positions内でamountが0の銘柄をnone_symsと定義したset型に入れる。
3.売りシグナルのループ文でnone_symsに該当する銘柄があったらcontinueする。

実装

サンプルコードは以下のリンクから参照してください。
https://factory.quantx.io/developer/82025c466c0e451ea4f6cd6d3cd26abb

また銘柄によってはwarningは出ますが出ている理由としては株を購入する際、売買条件にあっていないためです。

handle_signals関数以下を記述します。
def handle_signals(ctx, date, current):
  '''
  current: pd.DataFrame
  '''
  #市況シグナル
  market_sig = current["market:sig"]
  #利益確定及び損切りを格納するset型
  done_syms = set([])
  #portfolio.positionsに存在しない銘柄を格納するset型
  none_syms = set([])

  #ログ出力（確認用 日付、市況シグナル、portfolio.positions）
  # ctx.logger.debug(date)
  # ctx.logger.debug(market_sig)
  # ctx.logger.debug(ctx.portfolio.positions.items())

  #portfolio.positionsに銘柄が存在するかのチェック
  for (sym, val) in market_sig.items():
    judge = sym in ctx.portfolio.positions
    if (judge == False):
      none_syms.add(sym)
  for (sym, val) in ctx.portfolio.positions.items():
    c_amounts = val["amount"]
    if (c_amounts == 0):
      none_syms.add(sym)

  # 損切り、利確の設定    
  for (sym, val) in ctx.portfolio.positions.items():
    returns = val["returns"]
    amounts = val["amount"]
    if returns < -0.04:
      sec = ctx.getSecurity(sym)
      sec.order(-val["amount"], comment="損切り(%f)" % returns)
      done_syms.add(sym)
    elif returns > 0.05:
      sec = ctx.getSecurity(sym)
      sec.order(-val["amount"], comment="利益確定売(%f)" % returns)
      done_syms.add(sym)

  #ログ出力（確認用 保有していない銘柄、保有している銘柄）
  # ctx.logger.debug(none_syms)
  # ctx.logger.debug(done_syms)

  # 買いシグナル
  buy = market_sig[market_sig > 0.0]
  for (sym, val) in buy.items():
    if sym in done_syms:
      continue
    sec = ctx.getSecurity(sym)
    sec.order(sec.unit() * 1, orderType=maron.OrderType.MARKET_OPEN, comment="SIGNAL BUY")
    # ctx.logger.debug("購入: %s,  %f" % (sec.code(), val))
    pass

  # 売りシグナル
  sell = market_sig[market_sig < 0.0]
  for (sym, val) in sell.items():
    if sym in done_syms:
      continue
    if sym in none_syms:
      continue
    sec = ctx.getSecurity(sym)
    sec.order_target(0,orderType=maron.OrderType.MARKET_OPEN, comment="SIGNAL SELL")
    # ctx.logger.debug("売却: %s,  %f" % (sec.code(), val))
    pass
これで保有していない銘柄に対して売り注文を執行しないようになりました。

終わりに

portfolio_positions等の銘柄情報でかなり色々なデータを取ることができるとわかりました。
次回はこの辺のQuantX独自APIについてまとめようかなと思います。

宣伝

勉強会やってます！

日時：毎週金曜日19時〜
場所：神田　千代田共同ビル４階　SmartTrade社オフィス
内容：初心者(プログラミングってものを知らなくてもOK)向けに初心者(私とか)がこんな内容をハンズオン(一緒にやる事)で解説しています
備考：猛者の方も是非御鞭撻にいらして下さい、そして開発・伝導者になりましょう！

もくもく会もやってます！

日時：毎週水曜日18時〜
場所：神田　千代田共同ビル４階　SmartTrade社オフィス
内容：基本黙々と自習しながら猛者の方に質問して強くなっていく会
備考：お菓子と終わりにお酒を飲みながら参加者と歓談できます！

詳細はこちらです

週によっては開催されない週もあります。
また勉強会参加、もくもく会参加には基本的に事前に参加登録をしてください！
Pythonアルゴリズム勉強会HP：https://python-algo.connpass.com/
(connpassって言うイベントサイトに飛びます)

ストアもあります

システムトレードの開発者が作ったアルゴリズムがQuantX Storeで販売されています！
詳細は以下のリンクから
https://quantx.io/

免責注意事項

このコードや購入したアルゴリズム及び知識を使った実際の取引で生じた損益に関しては一切の責任を負いかねますので御了承下さい

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T14:41:41+09:00

マルチスレッドは良い

参考

いろんな記事でマルチスレッドを見てきたが、これが一番良い

https://news.mynavi.jp/article/python-33/

僕を含め、python初心者の方に知って欲しい

multi-threads

import threading, time, requests

def get_html(url):
    current_time = time.time()
    response = requests.get(url)
    html = response.text
    print(url + ': ' + str(time.time() - current_time))

urls = [
    'http://www.google.com',
    'http://www.yahoo.co.jp/',
    'https://www.bing.com/'
]
threads = []

# Start Threads
current_time = time.time()
for url in urls:
    thread = threading.Thread(target=get_html, args=(url,))
    thread.start()
    threads.append(thread)

# Wait Threads
for thread in threads:
    thread.join()

print('Time: ' + str(time.time() - current_time))

# Wait Threadsでスレッドの終了を待つ必要がある、コメントアウトして実行すればわかるが、最後のprintが先に実行されてしまう。

次に、継承によるマルチスレッド向けのクラスの作成と制御（lock）をメモする。
ちなみに、マルチスレッドの返り値取得方法をメモ

import threading, time

class MyThread(threading.Thread):
    def __init__(self, count):
        threading.Thread.__init__(self)
        self.count = count
        self.return_value = None   # RETURN VALUE

    def run(self):
        sum_value = 0
        for i in range(1, 1 + self.count):
            sum_value += i
            time.sleep(0.1)
        self.return_value = sum_value   # SET RETURN VALUE

    def get_value(self):  # GETTER METHOD
        return self.return_value

thread1 = MyThread(5)
thread1.start()
thread1.join()
print(thread1.return_value)  # 15
print(thread1.get_value())   # 15

個人的にはインスタンス変数を取得するメソッドを書く必要性がわからない。コメントください

次に、マルチスレッドの問題、コンフリクト的な問題、の例と解決策をメモ

商品Xの在庫が4あるとしましょう。店舗でお客さんがそれを購入する場合、実店舗のシステムAが在庫数を読み取り、その値からひとつ引いた3を新しい在庫数として登録します。ただ、そのシステムAが在庫数4を読み取り3を設定する間に、オンラインストアのシステムBでも同じ商品Xが購入されたとしましょう。オンラインストアも実店舗と同様に在庫数4を読み取り、新しく在庫数3を設定しようとします。このような状況では、上記のように実際には2つの製品が売れているにもかかわらず、在庫数は3となってしまいます。

import threading, time

resource = 5

class MyThread(threading.Thread):
    def __init__(self, name, sleep_time):
        threading.Thread.__init__(self)
        self.name = name
        self.sleep_time = sleep_time

    def run(self):
        global resource

        # read
        count = resource
        print(self.name + 'から見た現在の在庫は ' + str(resource))

        # do something
        print(self.name + 'が購入手続きなどを行っている')
        time.sleep(self.sleep_time)

        # write
        resource = count - 1
        print(self.name + 'から見た残りの在庫は ' + str(resource))


thread1 = MyThread('Aくん', 5)
thread2 = MyThread('Bくん', 3)
thread1.start()
time.sleep(1)
thread2.start()
thread1.join()
thread2.join()
print('最終的な在庫は ' + str(resource))

以下が出力です

Aくんから見た現在の在庫は 5
Aくんが購入手続きなどを行っている
Bくんから見た現在の在庫は 5
Bくんが購入手続きなどを行っている
Bくんから見た残りの在庫は 4
Aくんから見た残りの在庫は 4
最終的な在庫は 4

オンラインショッピングで上のようなことが起こると非常にまずい、ということで

threading.Lock() を使ってこの問題を解決する

あるスレッドがlockをしている場合は、別のスレッドは「lockする箇所」で待機し、リソースがreleaseされると、次は自分がlockをします。

import threading, time

resource = 5
global_lock = threading.Lock()   # LOCK OBJECT


class MyThread(threading.Thread):
    def __init__(self, name, sleep_time):
        threading.Thread.__init__(self)
        self.name = name
        self.sleep_time = sleep_time

    def run(self):
        global resource
        global global_lock

        # LOCK
        global_lock.acquire()

        # read
        count = resource
        print(self.name + 'から見た現在の在庫は ' + str(resource))

        # do something
        print(self.name + 'が購入手続きなどを行っている')
        time.sleep(self.sleep_time)

        # write
        resource = count - 1
        print(self.name + 'から見た残りの在庫は ' + str(resource))

        # RELEASE
        global_lock.release()


thread1 = MyThread('Aくん', 5)
thread2 = MyThread('Bくん', 3)
thread1.start()
time.sleep(1)
thread2.start()

thread1.join()
thread2.join()
print('最終的な在庫は ' + str(resource))

なお、出力

Aくんから見た現在の在庫は 5
Aくんが購入手続きなどを行っている
Aくんから見た残りの在庫は 4
Bくんから見た現在の在庫は 4
Bくんが購入手続きなどを行っている
Bくんから見た残りの在庫は 3
最終的な在庫は 3

これでスレッドAがlockしてからreleaseするまで、スレッドBはread、write処理ができない
また、

threading.Lockのインスタンス については
同じクラスのスレッドだけで同期をとるのであればクラス変数を利用することを推奨 らしい

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T14:37:26+09:00

【Django】テキストファイルを基にメーラー起動リンク作成

mailtoを使えば、a href～で簡単にリンク作れるんですけども、
リンクにカーソルを合わせたとき左下に～の部分が全部出てしまうんですよね。
メールアドレスだけならまだしも、本文がずらーっと出るのはダサいので嫌です！

よって、考えた手段がこちら
①txtから文章を読み込んで（python）
②読み込んだ文章とmailto使ってメーラー起動リンクを作成（javascript）

Djangoでこれを実現するには画面側からpython側にイベントを送らないといけないと思ったので、
javascriptやAjaxを調べつつチャレンジ。

javascript

こちらを参考に、クリックした箇所に対してメーラー起動リンクを作成するようなjsを書きます。
javascriptとはコンソールログで遊んだ程度の付き合いしかありません。

Django側からデータ受け取るためにAjax使ってますが、
こちらからは特に何も送ってないので変な使い方してる気はします。
function sendmail() {
    var is = this;

    $.ajax({
        url: "read_mail_text/",//Djangoのurlsに書いてあるパス
        type : 'POST',
        dataType: 'json',
        data : {app: null},//空っぽのjsonを送信
        cache: false,
        // 成功
        success: function(data) {
            //中身確認
            console.log(data["ac"]);
            //受け取った値からmailtoでメーラー起動用のリンク作成
            is.location.href = "mailto:" + data["ac"] + "?subject=" + data["sub"] + "&body=" + data["text"];

        },
        // エラー処理
        error: function(XMLHttpRequest, textStatus, errorThrown) {
            alert("Error : " + errorThrown);
       }
   });

}
Django

やることは
①Ajaxから呼び出されて
②txtから文章を読み込んで
③jsonを返す
……かと

①Ajaxにお呼ばれするために（Django）
path("read_mail_text/",[アプリ名].read_mail_text,name="read_mail_text")
こんな感じでパスを書いておきます。
敢えて書く必要もない気もしますが、一応メモ。

②txtの中身

最終的にmailtoに組み込むのでそれを準拠に、一行ずつ変数指定しました。
txtの中身は以下のような感じです
ac = "V3 <damidami@v3.mask>"; //送り先のアドレス
sub = "年賀状";//件名
text = "新年おめでとう。%0D%0A%0D%0A今年こそ、裏切り者は殺す。%0D%0AADKデストロンを代表して、ヨロイ元帥。";//本文
エスケープシーケンスのタグでも改行が出来なかった気がするので、改行コードがギエピーしてます。

③javascriptにjsonを返す（Django）

Ajaxのリクエストを受け取ってるので（たぶん）、引数にrequestを入れます。
現在地が分からなかったので、相対パス使ってます。
      def read_mail_text(self,request):
          #txt読込→変数に入れる
          path =os.path.dirname(__file__)+"\\mail_text.txt"

          myfile =File(open(path,  'rb'))
          txt=[]
          for t in myfile:
              txt.append(t.decode('utf-8'))

          #json形式で返す
          return JsonResponse({'ac':txt[0],'sub':txt[1],'text':txt[2]})
完成

自分の環境で動きましたのでひとまず、題名の内容はクリア。
汎用性が低いのがネックですね。

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T14:15:35+09:00

Pythonの関数について(乱数)

乱数について

乱数値（ランダムな値）を取得する方法

randomモジュールにはさまざまな乱数値の取得方法があります。用途に応じて適切なものを選ぶとよいでしょう。
import random

print(random.random())
print(random.uniform(1,100))
print(random.randint(1,100))
print(random.choice("1234567890"))

list=["a","b","c"]
random.shuffle(list)
print(list)
それぞれ異なる出力結果となります。

random.random()
0.0~1.0までのfloat値を取得します。

ramdom.uniform(x,y)
x~yまでのfloat値を取得します。

random.randint(x,y)
x~yまでのint値を取得します。

random.choice("1234567890")
1234567890の内から一つの要素をランダムで取り出します。

random.shuffle(list)
listの中身をシャッフルします。

【ホーム】

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T11:05:48+09:00

Arc GISに追加されているshapefileをQGISに（またはその逆）

Python（一部手作業）で、ArcMap上に追加されているレイヤをQGISに追加する（またはその逆の）手順を、紹介します。

Ⅰ.　ArcGISからQGISへ

Ⅰ-1. ArcGISのPythonコンソールでの作業

ArcMapのPythonコンソールに以下を実行します。
#mxdの定義
mxd = arcpy.mapping.MapDocument("CURRENT")
#get layer list
lyr = arcpy.mapping.ListLayers(mxd)


#レイヤの保存先パスをリストに格納
lyrPath =[]
for l in lyr:
    if l.isFeatureLayer ==True or l.isRasterLayer ==True:
        lyrPath.append(l.dataSource)
print(lyrPath)
実行結果として、追加されているレイヤの保存先パスがリストで取得される。
通常「arcpy.mapping.ListLayers(mxd)」を実行するとレイヤの一覧を取得できるが、下の図のようにレイヤーがグループになっていると、そのグループもリストに加わってしまい、次の処理でそれぞれのレイヤのパスを取得する段階でエラーになってしまう。
そこで、それぞれのレイヤーのタイプを調べて、ベクターレイヤもしくはラスターレイヤだったらlyrPathという別のリストに入れてあげるということをしています。

（ここでは、基盤地図情報から用意したshpと10mDEMがリスト化されています。）

下の画像の赤枠で囲った、実行結果の保存パスが表示されている部分（「[」から、「]」まで）を範囲選択してコピーしてください。

Ⅰ-2. QGISのPythonコンソールでの作業

今度はQGISのpythonコンソールを立ち上げて、以下のコードの「lyrPath = 」に続いて、先ほどArcGISのPythonコンソールでコピーした内容を貼り付けて実行する。
lyrPath =  コピーした内容を貼り付け
（今回の例の場合だと、こうなる）

続いて、以下を実行する。
lyrPath.reverse()
for i in lyrPath:
    if i[len(i)-4:] ==".shp":
        iface.addVectorLayer(i,'','ogr')
    else:
        iface.addRasterLayer(i,"")
完成！

iface.addVector(Raster)Layerでそれぞれのデータを追加していますが、そのままリスト順で追加すると、次々にレイヤのトップに追加されて、最終的に、ArcGISでのレイヤの並びとQGISのレイヤの並びが逆になってしまいます。

そこでlyrPath.reverse()でリストの並びを逆にしてあげています。
次のforループではパスの最後4文字が「.shp」ならベクターで追加、そうでないならラスターとして追加という感じです。
（.shpじゃないベクターデータならどうするかということはとりあえず無視しています）

Ⅱ. QGISからArcGISへ

Ⅱ-1.　QGISのPythonコンソールでの作業

QGISのPythonコンソールに以下を実行。
lyrlist = QgsProject.instance().mapLayers()

lyrPath =[]
for i in lyrlist:
    lyrPath.append(lyrlist[i].dataProvider().dataSourceUri())

print(lyrPath)
実行結果として、レイヤーのパスをリストで取得します。
QGISではレイヤ内にグループがあっても、リストには無視してくれます。
しかしこちらで気になるのが、パスを取得すると各レイヤのパスの後ろに「|layerid=0」というのがくっついています。（図中紫枠）
さらにフィルターをかけているレイヤには、後ろに「subset=フィルター式」がはいります（図中青枠）。
これらは、データをロードするときに邪魔になりますが、次の処理で対処するとしてとりあえず気にせず「print(lyrPath)」の実行結果部分（図中赤枠範囲）をコピーします。

Ⅱ-2.　ArcGISのPythonコンソールでの作業

今度はArcGISのpythonコンソールを立ち上げて、以下のコードの「lyrPath = 」に続いて、先ほどQGISのPythonコンソールでコピーしたレイヤのパスの内容を貼り付けて実行する。
lyrPath =  コピーした内容を貼り付け
（今回の場合だと、こうなる）

続いて以下を実行。
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd,"")[0]


for i in lyrPath:
    if i.find(".shp") !=-1:
        addlayer = arcpy.mapping.Layer(i[0:i.find(".shp") +4])
        arcpy.mapping.AddLayer(df, addlayer, "BOTTOM")
    else:
        addlayer = arcpy.mapping.Layer(i)
        arcpy.mapping.AddLayer(df, addlayer, "BOTTOM")
実行結果

無事追加される。

forループでlyrPath内をひとつずつ追加していくのですが、上で書いたとおり、QGISで取得した保存先パスの後ろに「layerid」とか「subset」とか余計な文字が入っています。
そこで、パス内に「.shp」という文字が入っているかを検索し、あればその.shpまでの文字数文までをレイヤーとして読み込む。
ラスターの時にはそういった余計な文字が入っていないので、「.shp」という文字が入っていなければ（＝ラスターの時は）、そのままレイヤとして読み込む。
ということです。
追加の際には"BOTTOM"を指定して、レイヤの最下層に追加するようにしてやれば、QGISでのレイヤの並びと同じになります。

Ⅲ. 課題

kmlやcsv、xyzタイルなどシェープファイルやラスター以外のデータがレイヤとしてあった場合を対応できてないです。

あと、QGISでレイヤのパスを取得するとフィルター情報も拾えるので、Arcに追加するときにそれも反映できればいいですね。

それはまた追々。

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T11:00:19+09:00

seabornで画像を保存する

バージョン

seaborn: 0.9.0

seabornで画像を保存する

seabornで画像を保存したい場合
line_plot = sns.lineplot()
figure = line_plot.get_figure()
figure.savefig(save_path)
のように、プロットオブジェクトから図オブジェクトを取り出し、そしてそれを保存すると良いです。

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T10:41:37+09:00

gensimのpoincare回すときの注意点（随時追加）

学習データは2つのアイテムのタプル
例：AとBというアイテムがあったとき，
train_date = [(A, B)]

model.saveしたときに，指定したファイルとhoge.model.kv.syn0.npyが出力されるが，このよくわからんファイルは指定したファイルと一緒において置かないとPoincareModel.loadできない

思いついたら随時更新します．

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T07:58:44+09:00

SimPyで離散事象シミュレーション(3) リソースを理解しよう：StoreとResource編

はじめに

SimPyというPythonの離散事象シミュレーション用のパッケージを見つけて試してみたら気に入ったので自分用の備忘録も兼ねて使い方をまとめていく．

SimPyのドキュメント

SimPyのソース

第1回の記事（はじめの一歩）

第2回の記事（リソースを理解しよう：Container編）

第3回目の今回は，「resources関連のモジュール群」の中で前回詳しく触れられなかったstore.pyとresource.pyの2つのモジュールをとりあげる．

Store系のリソース

store.pyモジュールの中にはStore系のリソース（Store，PriorityStore，FilterStoreの3つ）が定義されている．基本となるのがStoreクラスで，これはBaseResourceクラスを継承したサブクラスになっている．残る2つはStoreクラスを継承して少しカスタマイズしたものになっている．ここでは，基本となるStoreクラスを中心にみていこう．

Storeの概要

前回紹介したContainerクラスは，液体などの連続的なオブジェクト，もしくは均一で互いに区別できない離散的なオブジェクトを蓄える容器のようなリソースであった．そのため，容器の中身を単に数量（level）で把握しておけばよかった．

これに対して，Storeクラスは，離散的でかつ個別に区別できるアイテムを保持しておく容器のようなリソースである．したがって，容器の中身をアイテムのリスト（items）で把握している．それにあわせて，Putクラスを継承したStorePutクラスには，容器に格納したいアイテム（item）の情報が追加されている．

そして，_do_put()メソッドでは，Storeの容量に余裕があれば（itemsリストの長さがcapacity未満であれば），StorePut事象のitemをitemsリストの末尾に追加してから，ok=Trueで対応するStorePut事象をトリガーしている．

また，_do_get()メソッドでは，Store内にアイテムが存在すれば（itemsリストの長さが1以上であれば），itemsリストの先頭の要素を取り出し，それをvalueに代入し，ok=Trueで対応するStoreGet事象をトリガーしている．

これらの実装の結果，Storeクラスは，アイテムの先入れ先出し（FIFO）バッファ，あるいはQueueになっていることがわかる．仮に，これを後入れ先出し（LIFO）バッファ，あるいはStackに変更したいとすると，_do_get()メソッドをオーバライドして，リストの取り出し順を末尾からに変更すればいいだろう（ただし，有限のcapacityを設定してしまうと全体としては必ずしもStackにはならない）．

Storeクラスのインスタンスは，
simpy.resources.store.Store(env, capacity=float('inf'))
simpy.Store(env, capacity=float('inf'))
のいずれかで生成できる（2つ目はショートカット）．生成時に特に指定しないと，capacityは無限になる．

Storeの使用例

簡単な例として，makerプロセスが100個のproductを1つずつ生産して順にStoreに格納していき，buyerプロセスはそれらをStoreから1つずつ購入していく，という状況をモデル化してみた．
import random
import simpy

def maker(env):
    for i in range(100):
        time_to = random.expovariate(1)
        yield env.timeout(time_to)
        yield env.store.put('product_' +str(i))  # put string "product_i" into Store

def buyer(env):
    for j in range(100):
        time_to = random.expovariate(1)
        yield env.timeout(time_to)
        item = yield env.store.get()  # get string "product_i" from Store and assign it to item
        print('buyer_{} bought {}.'.format(j, item))

def main():
    env = simpy.Environment()
    env.store = simpy.Store(env)
    env.process(maker(env))
    env.process(buyer(env))
    env.run()

if __name__ == "__main__":
    main()
これを動かしてみると，StoreがFIFOになっていることが確認できる．次に，Storeを少しカスタマイズしてみよう．
import random
import simpy

class CustomStore(simpy.Store):
    def __init__(self, env, capacity=float('inf')):
        super(CustomStore, self).__init__(env, capacity)

    def _do_get(self, event):
        if self.items:
            event.succeed(self.items.pop(len(self.items) -1))

def maker(env):
    for i in range(100):
        time_to = random.expovariate(1)
        yield env.timeout(time_to)
        yield env.store.put('product_' +str(i))  # put string "product_i" into Store

def buyer(env):
    for j in range(100):
        time_to = random.expovariate(1)
        yield env.timeout(time_to)
        item = yield env.store.get()  # get string "product_i" from Store and assign it to item
        print('buyer_{} bought {}.'.format(j, item))

def main():
    env = simpy.Environment()
    env.store = CustomStore(env)
    env.process(maker(env))
    env.process(buyer(env))
    env.run()

if __name__ == "__main__":
    main()
最初に，Storeクラスを継承したCustomStoreというサブクラスを作成している．これは，itemsリストの末尾から順に要素を取得するように_do_get()メソッドをオーバーライドしたものである（FIFOをLIFOに変更）．

main()関数では，Storeの代わりにCustomStoreを利用して上とまったく同じ処理を行っている．これを実際に動かしてみると，先ほどとは購買されていく製品の順序が異なることが確認できる．

PriorityStoreとFilterStore

PriorityStoreは，itemsをリストではなくヒープにすることで，アイテムをなんらかの優先順位の順に取り出せるように，Storeクラスを拡張したものである．したがって，PriorityStoreに格納するアイテムはorderableでないといけない（__lt__()メソッドをもつこと）．そのために，PriorityItemクラスが用意されているので，必要に応じて利用しよう．

FilterStoreは，StoreGetクラスをさらに継承してFilterStoreGetクラスを定義し，それをGet事象として利用している．このFilterStoreGetクラスには，itemをフィルタリングするメソッド（filter()）が追加されていて，_do_get()でitemsリストから要素を抽出する際に，filter()がTrueを返すという条件を加えている．

Resource系のリソース

resource.pyモジュールの中にはResource系のリソース（Resource，PriorityResource，PreemptiveResourceの3つ）が定義されている．基本となるのがResourceクラスで，これはBaseResourceクラスを継承したサブクラスになっている．PriorityResourceはResourceを，PreemptiveResourceはPriorityResourceをそれぞれ継承して少しカスタマイズしたものになっている．ここでは，基本となるResourceクラスを中心にみていこう．

Resourceは，プロセスが使用する道具，機械などを表していて，同時に利用可能な個数に上限（capacity）があるため，取り合いが起こる．これに対するPutとGetの事象としては，Put/Getクラスを継承したサブクラスのRequest/Releaseクラスを用いるようになっている．

Resourceクラスのput_queueはqueueに改名されており，そのリストにRequest事象が追加されていく．リソースの中身はusersという名称のリスト（Storeクラスのitemsリストに対応）で表され，その中にもRequest事象が格納される．このusersリストの長さを返す属性countも追加されている．

解釈としては，usersリストに含まれるRequest事象（を送出したプロセス）が現在そのリソースを使用中であり，queueリストに含まれるRequest事象（を送出したプロセス）はリソースが入手可能になるのを待機中ということになる．

Releaseクラスは対応するRequest事象への参照をもち，インスタンス（Release事象）が生成されるとすぐにリソースの_do_get()メソッドが呼ばれる．そして，_do_get()メソッドは，それに対応するRequest事象をusersリストから削除して，このRelease事象をok=Trueでトリガーする．

Requestクラスは，yieldする際にwithでコンテキスト化しておくと自動的にrelease()を呼ぶ（対応するRelease事象を生成する）仕組みが追加されている（ので，通常はそのスタイルでコーディングすること）．これでリソースの使用中に依頼元のプロセスから取り消しが入った場合などにも自動的にrelease()が走るようになる（前回紹介したcancel()とrelease()の違いに注意する）．

Requestクラスの_do_put()メソッドは，usersリストの長さがcapacity未満であれば，Request事象をusersリストの末尾に追加し，その事象の（リソース使用開始時刻を表す）属性usage_sinceにenv.nowの値を代入する．

Resourceクラスのインスタンスは，
simpy.resources.resource.Resource(env, capacity=1)
simpy.Resource(env, capacity=1)
のいずれかで生成できる（2つ目はショートカット）．生成時に特に指定しないと，capacityは1になる．

Resourceの使用例

ここではよくあるジョブショップスケジューリングの例を挙げておこう．
import random
import simpy

MACHINE = 3  # number of machines
JOB = 5  # number of jobs
OPR_MIN = 2  # minimum number of operations
OPR_MAX = 5  # maximum number of operations

def job(env, j):
    n_o = random.randint(OPR_MIN, OPR_MAX)  # number of operations of job j
    p_r = random.choices(range(MACHINE), k=n_o)  # processing route
    p_t = [random.randint(2,7) for o in range(n_o)]  # processing times
    for o in range(n_o):
        m = env.machines[p_r[o]]  # which machine to use
        with m.request() as req:
            yield req
            s_t = round(env.now)  # start time
            yield env.timeout(p_t[o])
            c_t = round(env.now)  # completion time
            print('o_{}_{} on m_{}: {}-{}'.format(j, o, p_r[o], s_t, c_t))

def main():
    env = simpy.Environment()
    env.machines = [simpy.Resource(env) for m in range(MACHINE)]
    for j in range(JOB):
        env.process(job(env, j))
    env.run()

if __name__ == "__main__":
    main()
プロセスjob()の中で，まずそのジョブのオペレーション数と加工経路，各オペレーションの処理時間をランダムに設定している．そして，加工経路の順に機械をリクエストし，確保できたら専有し，加工時間が経過したらリリースするというサイクルを繰り返している．`

maim()関数の中を見ればわかるように，各機械はResourceでモデル化されている．したがって，上記のコードでは，すべての機械がディスパッチングルールとしてFIFOを使用していることになる．なお，下で紹介するPriorityResourceを使えば，これをSPTなどの他のディスパッチングルールに簡単に変更してみることもできる（ので，ぜひ試してみてほしい）．

PriorityResourceとPreemptiveResource

PriorityResourceは，queueを，なんらかのkeyでソート可能なSortedQueueに変更している．また，PriorityRequestというRequestのサブクラスを定義して，それをPut事象として用いている．これらによって，リソースを利用するプロセスの優先順位付け（優先度の高いプロセスから順にリソースが割り当てられていく）を可能にしている点が標準のResourceとの違いである．

PreemptiveResourceは，リソースを割り当てる際の優先順位付けに加えて，さらに現在リソースを使用中のプロセスを中断させて，優先度の高いプロセスが優先度の低いプロセスから使用中のリソースを奪い取る行為（Preempt）を考慮できるようにしている．このとき，PreemptされたPriorityRequest事象はusersリストから削除され，それに対応するプロセスにはInterrupt例外が投げ込まれる（ので，プロセス側で必要に応じてこの例外を処理する）．

まとめ

今回はSimPyの「resources関連のモジュール群」のうち，store.pyとresource.pyの内容についてまとめた．StoreクラスとResourceクラスの標準的な使い方も紹介したので，機会があればぜひ試してみてほしい．次回はシミュレーションの途中経過のアニメーション化についてまとめたいと思う．

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T07:43:48+09:00

Pythonチュートリアルではえ＾〜ってo思ったことまとめ

はじめに

現在Python公式チュートリアルで、Pythonの勉強を再び一からやっています。

そこで、自分がはえ＾〜って思ったことを随時まとめていこうと思います。

8章

例外の発生

Pythonでは、構文エラー自体は実行時前にエラーが出てわかる（そりゃまあ）

ただ、C++などと違って、コンパイルがなく、型のミスや明らかに間違えている部分を検出することができず、実行時にしか分からない。

そのため、
a = 0

print(20 / a)
Traceback (most recent call last):
  File "/Users/students/Desktop/研究実装/Python_Tutorial/main.py", line 5, in <module>
    print(20 / a)
ZeroDivisionError: division by zero
上記のようなZeroDivisionErrorを発動し、その場で実行が終了してしまう。

このような例外にはtry-except文を使用する。

例外の補足
while True:
    try:
        x = int(input())
        y = 30 / x

        print("here")

    except (ValueError, ZeroDivisionError, NameError) as e:
        print(e)
例外のキャッチには、Pythonではtry-except構文を利用する。
tryの中でエラーが発生すると、まずexceptに例外が投げられる。
ここで、一致するエラーがあれば、except内を実行し、try-exceptのあとに移行する。
そのため、上記の例では、ZeroDivisionErrorが発生すると、print('here')は実行されない。

また、もしexceptで補足されない場合さらに上の階層へエラーが投げられる。
そこでも補足されない場合、プログラムが停止する。
def f():
    try:
        x = int(input())
        y = 30 / x

        print("here")

    except (ValueError, NameError) as e:
        print(e)


while True:
    try:
        f()
    except ZeroDivisionError:
        print("ZeroDivision")
上記で0除算が発生すると、f関数内で補足されないのでよりインデントが上の上位で補足される。

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T07:43:48+09:00

Pythonチュートリアルではえ＾〜って思ったことまとめ

はじめに

現在Python公式チュートリアルで、Pythonの勉強を再び一からやっています。

そこで、自分がはえ＾〜って思ったことを随時まとめていこうと思います。

8章　例外

例外の発生

Pythonでは、構文エラー自体は実行時前にエラーが出てわかる（そりゃまあ）

ただ、C++などと違って、コンパイルがなく、型のミスや明らかに間違えている部分を検出することができず、実行時にしか分からない。

そのため、
a = 0

print(20 / a)
Traceback (most recent call last):
  File "/Users/students/Desktop/研究実装/Python_Tutorial/main.py", line 5, in <module>
    print(20 / a)
ZeroDivisionError: division by zero
上記のようなZeroDivisionErrorを発動し、その場で実行が終了してしまう。

このような例外にはtry-except文を使用する。

例外の補足
while True:
    try:
        x = int(input())
        y = 30 / x

        print("here")

    except (ValueError, ZeroDivisionError, NameError) as e:
        print(e)
例外のキャッチには、Pythonではtry-except構文を利用する。
tryの中でエラーが発生すると、まずexceptに例外が投げられる。
ここで、一致するエラーがあれば、except内を実行し、try-exceptのあとに移行する。
そのため、上記の例では、ZeroDivisionErrorが発生すると、print('here')は実行されない。

また、もしexceptで補足されない場合さらに上の階層へエラーが投げられる。
そこでも補足されない場合、プログラムが停止する。
def f():
    try:
        x = int(input())
        y = 30 / x

        print("here")

    except (ValueError, NameError) as e:
        print(e)


while True:
    try:
        f()
    except ZeroDivisionError:
        print("ZeroDivision")
上記で0除算が発生すると、f関数内で補足されないのでよりインデントが上の上位で補足される。

例外の連続と継承

例外は、exceptを並べることで複数用意することができる。
ただし、上から見ていって一個でもエラーと種類があっていたらそこで終わる。

また、エラーを受け取る
つまりexcept SomeErrorとかける、SomeErrorcはExceptionというPythonの組み込みクラスを継承したクラスである。
そして、もし発生したエラーと基底クラスが同じならその時点でエラーが発生する。

例えば、以下の場合
class B(Exception):
    pass

class C(B):
    pass

class D(C):
    pass

for cls in [B, C, D]:
    try:
        raise cls()
    except D:
        print("D")
    except C:
        print("C")
    except B:
        print("B")

for cls in [B, C, D]:
    try:
        raise cls()
    except B:
        print("B")
    except C:
        print("C")
    except D:
        print("D")
DCBBCDと出力される。
最後の３つがBCDとなるのは、Dの親の親がBであり、Dなどの子供のエラーが発生しても親であるBがそのエラーを補足するためである。

例外のワイルドカード

すべてを吸収できる例外構文がある。
ただ、これを使うとなんでもかんでも補足してしまうのでやめよう。
ただ、関数の最後の用意しておいて、親に補足してもらうために利用するのはありっぽい。
def f():
    try:
        x = int(input())
        y = 30 / x

        print("here")

    except (ValueError, NameError) as e:
        print(e)

# 投げますねぇ！
    except:
        raise


while True:
    try:
        f()
    except ZeroDivisionError:
        print("ZeroDivision")
例外のクリーンナップ

try-exceptの最後にfinallyといういつもの構文を使える。
これによって、必ず最後に行いたい処理（ファイルのクローズなど）を行うことができて嬉しい

9章　クラス　

名前空間

Pythonでは、変数などの名前とオブジェクトのメモリを対応させるのを名前空間と読んでいる。

この名前空間は、それぞれ入れ子構造になっていて
例えば

mathモジュール内で使用されているすべての関数とそのメモリの名前空間

math.gcd内で使用されている変数とその名前空間

ユーザ定義モジュールの変数とその名前空間

モジュール内のクラスの名前空間

のように入れ子構造になっている。

そして、これらの名前空間は互いに参照しないようになっているため、名前衝突の危険性を減らすことができる。

スコープ

スコープや変数などの名前が探しにいける範囲。
これは名前空間で閉じられる？
ただ、モジュールA・モジュールBがあったとき、これらはまたげないので、イメージ的には完全に分離されている。

名前空間の中にスコープがある、というイメージで認識している。

globalとnonlocal
def scope_test():

    def do_local():
        spam = "local spam"

    def do_nonlocal():
        nonlocal spam
        spam = "nonlocal spam"

    def do_global():
        global spam
        spam = "global spam"

    spam = "test"
    do_local()
    print(spam)

    do_nonlocal()
    print(spam)

    do_global()

scope_test()
print(spam)

'''
test
nonlocal spam
global spam
'''
pythonにはglobalとnonlocalキーワードがある。
関数などそれぞれの名前空間で通常は、一番小さい変数・スコープを使用する。
ただし、もしその変数名が無いなら親名前空間に同じ名前の変数がないか探しに行く。

ここで、もしグローバル宣言されている変数と同じ名前を関数で使用すると、一番小さいスコープがデフォルトのため、ローカル変数ができて、グローバル宣言されている変数には影響が出ない。
そこで、globalキーワードをつけるとグローバル宣言されてる変数を直接参照できる。

また、同様にglobalじゃないけど一個上の階層の変数がほしいなど、ローカル変数の宣言じゃないんです〜ってときは、nonlocalキーワードを使用しよう。

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T02:13:30+09:00

鉄は熱いうちにたたけ～pyATS使ってみた第2弾～

RobotFrameWorkで簡単な判定をさせてみた

モチベが高い今こそ自分を叩いて学習しようという安易な気持ちでタイトルをつけました
タイトルからは想像が一切つかない内容ですみません、
今回は短めの記事になります。

前回のおさらい

持てる力をすべて使いpyats動かしてみた
詳細は以下
https://qiita.com/snamon/items/78eb3e314bf8e5787719

構成/対象機器

サーバは前回と同じくCentos8になります。
インストール手順などは前回の記事を参考ください。

今回のSnadBoxこちら

予約無しでいつでも使うことができます（便利）

今回の対象機器はIOSXRになります
Cisco IOS XR Software, Version 6.5.3
Testbedの作成

今回はpyatsのprojectは作成せずにRobotファイルを作成して、showコマンドの判定をしてみます。
Testbedは今回も必要になるので作成していきます。
[root@localhost ~]# source pyATS/bin/activate
(pyATS) [root@localhost ~]#
(pyATS) [root@localhost ~]# mkdir Robot
(pyATS) [root@localhost ~]#
(pyATS) [root@localhost ~]# cd Robot/
(pyATS) [root@localhost Robot]# vi testbed.yml
そしてTestbed作成ですが、
その前に一度対象機器にTeratermなどでLoginしておきます
SandBoxに入りUser/passとportを確認してSSH
このときプロンプトに表示されるHostnameも確認しておきましょう

では今度こそ作成
testbed.yml
---
testbed:
  name: RobotFrameWork_Testbed
  tacacs:
      username: admin
  passwords:
      tacacs: C1sco12345

devices:
  cisco.test:
    connections:
      defaults:
        class: 'unicon.Unicon'
      vty:
        protocol: ssh
        ip: sbx-iosxr-mgmt.cisco.com
        port: 8181
    os: iosxr
    type: iosxr
ここまで来たらtestbedが正しいのか接続確認を実施します。
方法はいくつかあるようですが、今回はgenie shellコマンドでConnectionCheckします。
(pyATS) [root@localhost Robot]# genie shell --testbed testbed.yml
Welcome to Genie Interactive Shell
==================================
Python 3.6.8 (default, Oct  7 2019, 17:58:22)
[GCC 8.2.1 20180905 (Red Hat 8.2.1-3)]

>>> from genie.testbed import load
>>> testbed = load('testbed.yml')
------------------------------------------------------------------------------- 
>>> test = testbed.devices['cisco.test']
>>> test.connect()
[2019-10-23 12:06:17,190] +++ cisco.test logfile /tmp/cisco.test-cli-20191023T120617188.log +++
[2019-10-23 12:06:17,192] +++ Unicon plugin iosxr +++
The authenticity of host '[64.103.37.3]:8181 ([64.103.37.3]:8181)' can't be established.
RSA key fingerprint is SHA256:11wzj9NmCIPoZpMxGcUtYYYkvTozWnFnIF8eKaMIWoI.
Are you sure you want to continue connecting (yes/no)?
[2019-10-23 12:06:18,714] +++ connection to spawn: ssh -l admin 64.103.37.3 -p 8181, id: 139739676523936 +++
[2019-10-23 12:06:18,715] connection to cisco.test
yes
Warning: Permanently added '[64.103.37.3]:8181' (RSA) to the list of known hosts.
Password:


RP/0/RP0/CPU0:cisco.test#
[2019-10-23 12:06:21,778] +++ initializing handle +++
[2019-10-23 12:06:21,782] +++ cisco.test: executing command 'terminal length 0' +++
terminal length 0
Wed Oct 23 16:10:35.058 UTC
RP/0/RP0/CPU0:cisco.test#
[2019-10-23 12:06:22,424] +++ cisco.test: executing command 'terminal width 0' +++
terminal width 0
Wed Oct 23 16:10:35.683 UTC
RP/0/RP0/CPU0:cisco.test#
[2019-10-23 12:06:23,052] +++ cisco.test: config +++
configure terminal
Wed Oct 23 16:10:36.320 UTC
RP/0/RP0/CPU0:cisco.test(config)#no logging console
RP/0/RP0/CPU0:cisco.test(config)#line console
RP/0/RP0/CPU0:cisco.test(config-line)#exec-timeout 0 0
RP/0/RP0/CPU0:cisco.test(config-line)#absolute-timeout 0
RP/0/RP0/CPU0:cisco.test(config-line)#session-timeout 0
RP/0/RP0/CPU0:cisco.test(config-line)#line default
RP/0/RP0/CPU0:cisco.test(config-line)#exec-timeout 0 0
RP/0/RP0/CPU0:cisco.test(config-line)#absolute-timeout 0
RP/0/RP0/CPU0:cisco.test(config-line)#session-timeout 0
RP/0/RP0/CPU0:cisco.test(config-line)#commit
Wed Oct 23 16:10:42.953 UTC
RP/0/RP0/CPU0:cisco.test(config-line)#end
RP/0/RP0/CPU0:cisco.test#
"Are you sure you want to continue connecting (yes/no)? yes\r\nWarning: Permanently added '[64.103.37.3]:8181' (RSA) to the list of known hosts.\r\r\nPassword: \r\n\r\n\r\nRP/0/RP0/CPU0:cisco.test#"
>>>
問題なさそうですね。

今回のテスト内容

RobotFrameWorkを使って以下をチェックします。

10.10.20.0/24のルーティングテーブルに登録されていること

デフォルトルートがルーティングテーブルに登録されていること

Gigabit0/0/0/0がAdminDownであること

さっそくテスト用のファイルを作成します
NWtest.robot
**** Settings ***
Library        ats.robot.pyATSRobot
Library        genie.libs.robot.GenieRobot
Library        unicon.robot.UniconRobot
Suite setup    Setup

*** Variables ***
${testbed}    testbed.yml
${ip1}  10.10.20.0/24

*** Test Cases ***
Verify Route
    ${cmd}=    parse "show route ipv4" on device "cisco.test"
    should contain  ${cmd}[vrf][default][address_family][ipv4][routes]  ${ip1}
    should contain  ${cmd}[vrf][default][address_family][ipv4][routes]  0.0.0.0/0

Verify Interface
    ${cmd}=    parse "show interface brief" on device "cisco.test"
    should contain  ${cmd}[interface][ethernet][GigabitEthernet0/0/0/0][status]  admin-down

*** Keywords ***
Setup
    use genie testbed "${testbed}"
パッと見た感じで多少無理がある感がでているがご愛嬌
コマンドは下記のgenieのドキュメントから選んできました。
https://developer.cisco.com/docs/genie-docs/

なーんかイケてない書き方なきがしてならないです。。。
実行ログ
------------------------------------------------------------------------------
NWtest                                                                | FAIL |
Suite setup failed:
ConnectionError: failed to connect to cisco.test
Failed while bringing device to "any" state

2 critical tests, 0 passed, 2 failed
2 tests total, 0 passed, 2 failed
==============================================================================
Output:  /root/Robot/output.xml
Log:     /root/Robot/log.html
Report:  /root/Robot/report.html
(pyATS) [root@localhost Robot]#
あれれ、、失敗した。。
ConnectionErrorですね、、何度か接続確認はできたので妙だなと思いつつも実機へ再度SSHしてみると。。。

ええぇ、、hostname変わってる。。
こういうところも自動化の注意点ですね、手で触るユーザがのこるとこんなリスクもあるのか。
今回はフリーの検証環境なのでしょうがないですね

hostnameを入れ替えて再度実行！
実行ログ
(pyATS) [root@localhost Robot]# robot NWtest.robot
==============================================================================
NWtest
==============================================================================
[ WARN ] Could not load the Datafile correctly
Verify Route                                                          | PASS |
------------------------------------------------------------------------------
Verify Interface                                                      | PASS |
------------------------------------------------------------------------------
NWtest                                                                | PASS |
2 critical tests, 2 passed, 0 failed
2 tests total, 2 passed, 0 failed
==============================================================================
Output:  /root/Robot/output.xml
Log:     /root/Robot/log.html
Report:  /root/Robot/report.html
(pyATS) [root@localhost Robot]#
無事成功！
実行結果として以下が作成されています（failのときも作成されてますね）

output.xml

log.html

report.html

このlog.htmlを表示すると.....

うん、いい感じですね
更に下の[+]を確認してきましょう

こんな感じです。
何やら確認をしてくれているようなんですが、、いまいちわからない。。

まとめ

うまくいってんだかよくわからねぇ
が、確認項目をミリミリかければ正常性確認はできそうな気がしました。
単体試験や結合試験でどうやって使っていくか？についてはこれからも考えていこうと思います。

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T01:36:57+09:00

FreePBX から slack 連携

FreePBXでslackに着信番号通知

これは

元職場で同僚だったmahさんのfax2slack をヒントに3CXにぶら下げた FreePBX で着信した電話番号の通知番号をslackに投下する奴です。

必要なもの

SIPサーバ(なんでも)

着信専用のFreePBX

Slack の WEB hook URL

FreePBXの設定(この手順わすれそうなのでこれを書いた)

WEBコンソールにログイン後アドミン → カスタム宛先を選択

Add Destination を押す

Target に mycustom-app,s,1 と入力し説明を適当に入れて Submit
/etc/asterisk/extensions_custom.conf
[mycustom-app]
exten => s,1,System("[SCRIPT PATH]" ${CALLERID(num)})
exten => s,n,Wait(300)
exten => s,n,Hangup
これで script path のスクリプトを実行してくれます。
user は asterisk の権限で動くので権限周りは確認しておくとよい。

あとはスクリプトを適当に書けば着信した番号を拾う事のできるシステムの完成。

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T01:36:54+09:00

[pandas]挙動が理解しがたいけど便利なgroupbyの使用例（複数カラムでの条件判定など）

はじめに

挙動を説明できるほど理解してないが、なんだかんだ便利なpandasのgroupbyの使用例をいくつか記載しておきます（備忘＆自戒を兼ねて）

サンプルデータ作成

結果

# check
df
"""
   x1  x2  gender young_old
0   3   5  female     young
1   8  10    male       old
2   2   7    male       old
3   9   0  female       old
4  11   6  female     young
5   1   4    male     young
"""

コード

import pandas as pd
import numpy as np

# 値
arr = np.arange(12)

# セグメント用
gender = ['male','female'] * 3
young_old = ['young','old'] * 3

# シャッフル
np.random.shuffle(arr)
np.random.shuffle(gender)
np.random.shuffle(young_old)

# 値のデータフレーム
df = pd.DataFrame(data=arr.reshape(6,2),
                  columns=['x1','x2']
                  )
# セグメントデータ追加
df['gender'] = gender
df['young_old'] = young_old

よく分からないが便利なgroupbyコード例

上から難しい順に3つ記載

よくわかんないけど便利なgroupby その3 複数カラム使って条件分岐

結果

# 他カラムからage_genderを生成した
# check
tmp[['gender','young_old','age_gender']]
"""
   gender young_old    age_gender
0  female     young  young_female
1    male       old      old_male
2    male       old      old_male
3  female       old    old_female
4  female     young  young_female
5    male     young    young_male
"""

コード

# これはかなり無茶なやり方・・・
def my_func_switch(xdf):
    # 1行になる時だけちゃんと動く これ外して動いても意図した計算にならないと思うのでやらない方がいい
    if xdf.shape[0]>1:raise('xdf.shape[0]==1のみ有効')
    # 条件分岐
    if xdf['gender'].values[0]=='male':
        if xdf['young_old'].values[0]=='young':
            return 'young_male'
        elif xdf['young_old'].values[0]=='old':
            return 'old_male'
    elif xdf['gender'].values[0]=='female':
        if xdf['young_old'].values[0]=='young':
            return 'young_female'
        elif xdf['young_old'].values[0]=='old':
            return 'old_female'

tmp = df.copy()  # なくてもいい。ここではdfを書き換えたくないのでやってる
tmp['dummy'] = np.arange(6)  # ユニークな値を持つカラムを作る
tmp['age_gender'] = (tmp
                       .groupby('dummy')
                       .apply(my_func_switch)
                       .values
                       )

よくわかんないけど便利なgroupby その2 データフレームを受け取って複数カラム使う

結果

# flgカラムを生成（x1<x2の判定結果）
# check
tmp[['x1', 'x2', 'flg']]
"""
   x1  x2    flg
0   3   5   True
1   8  10  False
2   2   7  False
3   9   0   True
4  11   6   True
5   1   4   True
"""

コード

tmp = df.copy()  # なくてもいい。ここではdfを書き換えたくないのでやってる
tmp['flg'] = (tmp
               .groupby('gender')
               .apply(lambda xdf: xdf['x1'] < xdf['x2'])
               .values
               )

よくわかんないけど便利なgroupby その1 x/x.max()みたいな

結果

# genderごとに最大値を算出し、最大値に対する割合のカラムを生成
# check
tmp[['x1','gender','genderごとのmaxに対する割合']]
"""
  x1  gender  genderごとのmaxに対する割合
0   3  female            0.272727
1   8    male            1.000000
2   2    male            0.250000
3   9  female            0.818182
4  11  female            1.000000
5   1    male            0.125000
"""

コード

a = (df
     .groupby('gender')
     ['x1']
     .apply(lambda x: x/ x.max())
     )
# check
a
"""
0    0.272727
1    1.000000
2    0.250000
3    0.818182
4    1.000000
5    0.125000
Name: x1, dtype: float64
"""

tmp = df.copy()
tmp['genderごとのmaxに対する割合'] = a

シンプルなgroupbyコード例

上から難しい順に6つ記載

# シンプル？なgroupby  その6 shift
tmp = df.copy()
tmp['shifted'] = (tmp
                 .groupby('gender')
                 ['x1']
                 .shift(-1)
                 .values
                 )
# check
tmp[['x1','gender','shifted']]
"""
   x1  gender  shifted
0   3  female      9.0
1   8    male      2.0
2   2    male      1.0
3   9  female     11.0
4  11  female      NaN
5   1    male      NaN
"""


# シンプルなgroupby　その5 apply&lambda利用
(df
 .groupby('gender')
 ['x1']
 .apply(lambda srs: srs.max())
 .reset_index()
)
# check
"""
   gender  x1
0  female  11
1    male   8
"""


# シンプルなgroupby　その4 apply&関数利用
def my_func_max(srs):
    # 型はSeries
    return srs.max()

(df
 .groupby('gender')
 ['x1']
 .apply(my_func_max)
 .reset_index()
)
# check
"""
   gender  x1
0  female  11
1    male   8
"""
# おまけ） 計算前にどのように値が分けられるか確認
# check
df.groupby('gender').groups
"""
{'female': Int64Index([0, 3, 4], dtype='int64'),
 'male': Int64Index([1, 2, 5], dtype='int64')}
"""


# シンプルなgroupby その3 agg
(df
 .groupby('gender')
 .agg({'x1':np.max, 'x2':np.max})
 .reset_index()
)
# check
"""
 gender  x1  x2
0  female  11   6
1    male   8  10
"""


# シンプルなgroupby その2
(df
 .groupby(['gender','young_old'])
 ['x1']
 .max()
 .reset_index()
)
# check
"""
   gender young_old  x1
0  female       old   9
1  female     young  11
2    male       old   8
3    male     young   1
"""


# シンプルなgroupby
(df
 .groupby('gender')
 ['x1']
 .max()
 .reset_index()
)
# check
"""
gender  x1
0  female  11
1    male   8
"""

参考ページ

groupbyの基本はこちらに良くまとまっています。素敵なページです。
https://qiita.com/propella/items/a9a32b878c77222630ae

終わり

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T01:36:54+09:00

[pandas]よくわかんないけど便利なgroupbyの応用例

はじめに

挙動が分かりづらいが、なんだかんだ便利な時があるpandasのgroupbyのコード例をいくつか記載しておきます（備忘兼ねて）

サンプルデータ作成

結果

# check
df
"""
   x1  x2  gender young_old
0   3   5  female     young
1   8  10    male       old
2   2   7    male       old
3   9   0  female       old
4  11   6  female     young
5   1   4    male     young
"""

コード

import pandas as pd
import numpy as np

# 値
arr = np.arange(12)

# セグメント用
gender = ['male','female'] * 3
young_old = ['young','old'] * 3

# シャッフル
np.random.shuffle(arr)
np.random.shuffle(gender)
np.random.shuffle(young_old)

# 値のデータフレーム
df = pd.DataFrame(data=arr.reshape(6,2),
                  columns=['x1','x2']
                  )
# セグメントデータ追加
df['gender'] = gender
df['young_old'] = young_old

よくわかんないけど便利なgroupbyコード例

上から難しい順に3つ記載

その3 複数カラム使って条件分岐

結果

# 他カラムからage_genderを生成した（単純な文字列結合処理ではない）
# check
tmp[['gender','young_old','age_gender']]
"""
   gender young_old    age_gender
0  female     young  young_female
1    male       old      old_male
2    male       old      old_male
3  female       old    old_female
4  female     young  young_female
5    male     young    young_male
"""

コード

def my_func_switch(xdf):
    # 1行になる時だけちゃんと動く これ外して動いても意図した計算にならないと思うのでやらない方がいい
    if xdf.shape[0]>1:raise('xdf.shape[0]==1のみ有効')
    # 条件分岐
    if xdf['gender'].values[0]=='male':
        if xdf['young_old'].values[0]=='young':
            return 'young_male'
        elif xdf['young_old'].values[0]=='old':
            return 'old_male'
    elif xdf['gender'].values[0]=='female':
        if xdf['young_old'].values[0]=='young':
            return 'young_female'
        elif xdf['young_old'].values[0]=='old':
            return 'old_female'

tmp = df.copy()  # なくてもいい。ここではdfを書き換えたくないのでやってる
tmp['dummy'] = np.arange(6)  # ユニークな値を持つカラムを作る
tmp['age_gender'] = (tmp
                       .groupby('dummy')
                       .apply(my_func_switch)
                       .values
                       )

その2 データフレームを受け取って複数カラム使う

結果

# flgカラムを生成（x1<x2の判定結果）
# check
tmp[['x1', 'x2', 'flg']]
"""
   x1  x2    flg
0   3   5   True
1   8  10  False
2   2   7  False
3   9   0   True
4  11   6   True
5   1   4   True
"""

コード

tmp = df.copy()  # なくてもいい。ここではdfを書き換えたくないのでやってる
tmp['flg'] = (tmp
               .groupby('gender')
               .apply(lambda xdf: xdf['x1'] < xdf['x2'])
               .values
               )

その1 x/x.max()みたいな

結果

# genderごとに最大値を算出し、最大値に対する割合のカラムを生成
# check
tmp[['x1','gender','genderごとのmaxに対する割合']]
"""
  x1  gender  genderごとのmaxに対する割合
0   3  female            0.272727
1   8    male            1.000000
2   2    male            0.250000
3   9  female            0.818182
4  11  female            1.000000
5   1    male            0.125000
"""

コード

a = (df
     .groupby('gender')
     ['x1']
     .apply(lambda x: x/ x.max())
     )
# check
a
"""
0    0.272727
1    1.000000
2    0.250000
3    0.818182
4    1.000000
5    0.125000
Name: x1, dtype: float64
"""

tmp = df.copy()
tmp['genderごとのmaxに対する割合'] = a

シンプルなgroupbyコード例

上から難しい順に6つ記載

# シンプル？なgroupby  その6 shift
tmp = df.copy()
tmp['shifted'] = (tmp
                 .groupby('gender')
                 ['x1']
                 .shift(-1)
                 .values
                 )
# check
tmp[['x1','gender','shifted']]
"""
   x1  gender  shifted
0   3  female      9.0
1   8    male      2.0
2   2    male      1.0
3   9  female     11.0
4  11  female      NaN
5   1    male      NaN
"""


# シンプルなgroupby　その5 apply&lambda利用
(df
 .groupby('gender')
 ['x1']
 .apply(lambda srs: srs.max())
 .reset_index()
)
# check
"""
   gender  x1
0  female  11
1    male   8
"""


# シンプルなgroupby　その4 apply&関数利用
def my_func_max(srs):
    # 型はSeries
    return srs.max()

(df
 .groupby('gender')
 ['x1']
 .apply(my_func_max)
 .reset_index()
)
# check
"""
   gender  x1
0  female  11
1    male   8
"""
# おまけ） 計算前にどのように値が分けられるか確認
# check
df.groupby('gender').groups
"""
{'female': Int64Index([0, 3, 4], dtype='int64'),
 'male': Int64Index([1, 2, 5], dtype='int64')}
"""


# シンプルなgroupby その3 agg
(df
 .groupby('gender')
 .agg({'x1':np.max, 'x2':np.max})
 .reset_index()
)
# check
"""
 gender  x1  x2
0  female  11   6
1    male   8  10
"""


# シンプルなgroupby その2
(df
 .groupby(['gender','young_old'])
 ['x1']
 .max()
 .reset_index()
)
# check
"""
   gender young_old  x1
0  female       old   9
1  female     young  11
2    male       old   8
3    male     young   1
"""


# シンプルなgroupby
(df
 .groupby('gender')
 ['x1']
 .max()
 .reset_index()
)
# check
"""
gender  x1
0  female  11
1    male   8
"""

参考ページ

groupbyの基本はこちらに良くまとまっています。素敵なページです。
https://qiita.com/propella/items/a9a32b878c77222630ae

終わり

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T01:19:38+09:00

Django フォームセットでPOSTでの取得の仕方がわからない。

form.py

class 案内Form(forms.Form):    
    会社別=forms.CharField(required=False,max_length=3,widget=forms.TextInput(attrs={'class': "form-control"}))

views.py

class 案内View(TemplateView):  
    def __init__(self):
        案内FormSet=formset_factory(案内Form,extra=6)
        #buf=フォーム初期値入力
        buf={'会社別':1,'会社別':2,'会社別':3,'会社別':4,･･･}
        formset =案内FormSet(initial=buf)
        self.params = {
            'guide':formset,                       
        }
    def post(self, request):           
        formset = 案内FormSet(request.POST)
        #print(formset)でエラーが出る。
        print(formset)        
        return render(request, 'jw/guide.html', self.params)

guide.html

 <tbody>
・・・・・省略・・・・ 
<tr>
  <td>1</td>
  <td><input type="text" name="form-0-会社別" value="1" class="form-control" maxlength="3" id="id_form-0-会社別"></td>
  <td>
    <ruby>
      <rb>十九</rb>
      <rt>ジュウク</rt>
      　<rb>太郎</rb>
      <rt>タロウ</rt>
    </ruby>
  </td>
  <td>10/04</td>
  <td>新規</td>
  <td>SN-1F</td>
  <td>提出物</td>
  <td>入金</td>
  <td></td>
  <td>株式会社ジューク</td>
</tr>
<tr>
  <td>2</td>
  <td><input type="text" name="form-1-会社別" value="2" class="form-control" maxlength="3" id="id_form-1-会社別"></td>
  <td>
    <ruby>
      <rb>山崎</rb>
      <rt>ヤマサキ</rt>
      　<rb>太郎</rb>
      <rt>タロウ</rt>
    </ruby>
  </td>
  <td>10/04</td>
  <td>新規</td>
  <td>SN-1F</td>
  <td>提出物</td>
  <td>入金</td>
  <td></td>
  <td>株式会社ヤマダ</td>
</tr>
・・・・・省略・・・・              
</tbody>

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T01:09:54+09:00

Survival Analysis with Python

備忘録としてまとめます。

生存分析は、生物の死や機械システムの故障など、あるイベントが発生するまでの期間を分析するための統計学的な手法の1つです。
ある時期を過ぎて生存する人の割合はどのくらい？あるイベントが起きた人のうち、その人たちはどのくらいの割合で天寿をまっとうする？ある特性は、生存の確率をどのように増加または減少させる？このような疑問に答える手法を模索するための学問ですね。

Kaplan-Meier法を用いた予測

ある時期に発生するイベントの割合から発生確率を予測します。
データはこちらを使います。
「Tongue」
https://drive.google.com/file/d/1fpiIaKFMc9QPbOamKdxO1sWROzf1DP1Z/view?usp=sharing
元データ

sksurvを使う
#sksurvを使う方法
!pip install scikit-survival
#してから
#CSVファイルをこのノートへロードします。
#ファイルを開くダイアログを表示
from google.colab import files
file = files.upload()

#ファイルをCSVとして読み込みます
import pandas as pd
import io
data = pd.read_csv(io.StringIO(file['tongue.csv'].decode('utf-8')),header=0)
#中身の大きさと先頭5行のデータの内容を確認してみます。
print('データの形式：{}'.format(data.shape))
print(data.head())

#イベントのデータ型がintだと駄目なようなので、真偽値に変換しておきます。
#一先ず、グループ１のみを計算対象にします。
f = data.type==1
T = data[f]['time']
event = data[f]['delta']
eventBool = []
for i in event:
  if(event[i]==1):
    eventBool.append(True)
  else :
    eventBool.append(False)

#予測します。
from sksurv.nonparametric import kaplan_meier_estimator
x, y = kaplan_meier_estimator(eventBool, time)

#グラフにします
import matplotlib.pyplot as plt
plt.step(x, y, where="post")
plt.ylim(0, 1)
plt.show()
lifelinesを使う
!pip install lifelines # このライブラリも便利です。https://lifelines.readthedocs.io/en/latest/Quickstart.html

#からの予測器を初期化して
from lifelines import KaplanMeierFitter
kmf = KaplanMeierFitter()

#真偽値変換無しでそのまま予測器の引数に渡します
f = data.type==1
T = data[f]['time']
C = data[f]['delta']

kmf.fit(T, event_observed=C)

#グラフにします
kmf.plot(title='Tumor DNA Profile 1')
比較

sksurvでも出来ますが、lifelinesライブラリを使って、グループ2（つまりtype:2）と比較をしてみます。
f2 = data.type==2
T2 = data[f2]['time']
C2 = data[f2]['delta']

ax = plt.subplot(111)

kmf.fit(T, event_observed=C, label=['Type 1 DNA'])
kmf.survival_function_.plot(ax=ax)
kmf.fit(T2, event_observed=C2, label=['Type 2 DNA'])
kmf.survival_function_.plot(ax=ax)

plt.title('Lifespans of different tumor DNA profile')
logrank-test

生存時間解析において、2群の生存時間に差があるかを検定するノンパラメトリックな検定手法のこと。
一般化Wilcoxon検定に対して相対的に後期に起きた死亡を重く評価するという特徴を持つ。時間経過と共に生存率曲線の差が開いてくるような場合、ログランク検定は一般化Wilcoxon検定に比べて検出力が高くなる。（https://bellcurve.jp/statistics/glossary/719.html）
from lifelines.statistics import logrank_test
summary_= logrank_test(T, T2, C, C2, alpha=99)
summary_.print_summary()
出力
t_0 = -1 null_distribution = chi squared degrees_of_freedom = 1 alpha = 99 test_statistic p -log2(p) 2.79 0.09 3.40

ハザード比の算出(Nelson-Aalen kernelを利用)

ハザード比も算出できます。
ある臨床試験で検討したい新薬Aと比較対象の薬剤Bとを比べたとき、ハザード比が1であれば2つの治療法に差はなく、ハザード比が1より小さい場合には治療Aの方が有効と判定され、その数値が小さいほど有効であるとされます。　例えばA薬と対象のB薬を比較するというある臨床試験でハザード比が0.94という結果であれば、A薬はB薬よりリスクを6%減少させたという意味になります。(https://oncolo.jp/dictionary/hr)
from lifelines import NelsonAalenFitter

naf = NelsonAalenFitter()
naf.fit(T, event_observed=C)

naf.plot(title='Nelson-Aalen Estimate')
最後に「Time Dependent AUC」

よくバイナリの予測器を作成する際にROC曲線を使いますが、このときのROCはあくまでもその時点で収集できたデータで作成した予測器の精度を示すROCです。
このROCを一定の期間ごとに計算し、それぞれの時点でデータを整えて分類モデルでROC（その時点の生存の正解を利用）を描き、AUCを計算し、経時的なAUCの変化を可視化したのがTime dependent AUCです。

これはsksurvライブラリを使って試してみます。
免疫グロブリン遊離L鎖κ/λ比のサンプルデータをロードします。
これは、以下のようなデータです。
7874サンプルで9つの項目があります。
age：年数
sex：F =女性、M =男性
sample.yr：血液サンプルが取得された暦年
kappa：FLCカッパ
lambda：FLCラムダ
flc.grp：被験者のFLCグループ
creatinine：血清クレアチニン
mgus：被験者がmonoclonal gammapothy（MGUS）と診断されたかどうか
chapter：国際疾病コードICD-9の章の見出しによる主な死因のグループ化
#データをロードします
from sksurv.datasets import load_flchain
x, y = load_flchain()

#今回はテスト用にデータを分けることにします
from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2, random_state=0)

#'age', 'creatinine', 'kappa', 'lambda'のみ学習に利用します
num_columns = ['age', 'creatinine', 'kappa', 'lambda']

#欠測を簡易的に単一補完します
from sklearn.impute import SimpleImputer
imputer = SimpleImputer().fit(x_train.loc[:, num_columns])
x_train = imputer.transform(x_train.loc[:, num_columns])
x_test = imputer.transform(x_test.loc[:, num_columns])

#正解ラベルからトレーニングとテストのそれぞれのデータセットの生存期間の最大最小を取得します。
y_events = y_train[y_train['death']]
train_min, train_max = y_events["futime"].min(), y_events["futime"].max()

y_events = y_test[y_test['death']]
test_min, test_max = y_events["futime"].min(), y_events["futime"].max()

#データの整合性を確認します。テストデータに用いるデータの生存期間はトレーニングに用いるデータの生存期間内になければなりません。
assert train_min <= test_min < test_max < train_max, \
    "time range or test data is not within time range of training data."

#今回は生存期間が、5から81までの間で15間隔で計算します。
import numpy as np
times = np.percentile(y["futime"], np.linspace(5, 81, 15))
print(times)

#time dependent aucを描きます。
from sksurv.metrics import cumulative_dynamic_auc
from sksurv.metrics import concordance_index_ipcw
import matplotlib.pyplot as plt

def plot_cumulative_dynamic_auc(risk_score, label, color=None):
    auc, mean_auc = cumulative_dynamic_auc(y_train, y_test, risk_score, times)

    plt.plot(times, auc, marker="o", color=color, label=label)
    plt.xlabel("days from enrollment")
    plt.ylabel("time-dependent AUC")
    plt.axhline(mean_auc, color=color, linestyle="--")
    plt.legend()

for i, col in enumerate(num_columns):
    plot_cumulative_dynamic_auc(x_test[:, i], col, color="C{}".format(i))
    ret = concordance_index_ipcw(y_train, y_test, x_test[:, i], tau=times[-1])
上記のグラフから、予測力高い特徴量は「Age, kappa/lambda, creatinine」の順であることがわかります。

以上です。

p.s

さらに、sksurvの開発者が深層学習を利用した生存分析のアプローチも紹介しています。
勉強しなくては。。
SurvivalAnalysisForDeepLearning

Reference

https://en.wikipedia.org/wiki/Survival_analysis

https://plot.ly/python/v3/ipython-notebooks/survival-analysis-r-vs-python/

https://k-d-w.org/slides/pyconuk-2017/

http://sinhrks.hatenablog.com/entry/2014/11/01/223504

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T00:53:20+09:00

Python3で深さ優先探索と幅優先探索の実装

幅優先探索と深さ優先探索をPythonでシンプルに実装してみました。
今回は以下のようなグラフを例にノードb（Start）からノードg（End）への探索を実装してみます。

入力

以下のようなDictionary形式で各ノードと隣接するノードを表すことにします。

search.py

graph = {
    'a': ['c'],
    'b': ['c','f'],
    'c': ['a','b','d'],
    'd': ['e','f'],
    'e': ['d'],
    'f': ['b','d','g'],
    'g': ['f']
}

深さ優先探索(Depth-First Search)

考え方

Start(b)の隣接ノードそれぞれについてEnd(g)と一致しているかチェック
【一致しない】→Startの隣接ノードをStartにして探索の関数を再帰する
【一致する】→終了

（mark_dictに既出のノードを記録して無限ループを回避している）

実装

search.py

import sys

# 深さ優先探索
def DFS(graph,start,end,mark_dict):
    mark_dict[start] = True
    for target in graph[start]:
        if target == end:
            print('Exist')
            sys.exit()
        elif target in mark_dict:
            pass
        else:
            mark_dict[target] = True
            DFS(graph,target,end,mark_dict)

DFS(graph,'b','g',{})

幅優先探索(Breadth-First Search)

考え方

Start(b)の隣接ノードがEnd(g)と一致しているかチェック
【一致しない】→Queueにその隣接ノードに隣接するノードを入れ、Queueの最初の要素をStartにしてループさせる
【一致する】→終了

実装

search.py

# queueの実装
class Queue:
    def __init__(self):
        self.queue = []
    # 先頭にデータを追加
    def add(self,x):
        self.queue.append(x)
    # 最後のデータを削除
    def remove(self):
        del self.queue[0]
    # 最後のデータを返す
    def peek(self):
        return self.queue[0]
    # 空かどうか確認
    def is_empty(self):
        if len(self.queue) == 0:
            return True
        else:
            return False

# 幅優先探索
def BFS(graph,start,end):
    mark_dict = {}
    queue = Queue()
    queue.add(start)
    while queue.is_empty() == False:
        target = queue.peek()
        queue.remove()
        if target in mark_dict:
            pass
        else:
            mark_dict[target] = True
            if target != end:
                for elem in graph[target]:
                    queue.add(elem)
            else:
                print('Exist')
                sys.exit()

BFS(graph,'b','g')

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T00:53:20+09:00

Python3で深さ優先探索と幅優先探索と双方向探索の実装

幅優先探索と深さ優先探索をPythonでシンプルに実装してみました。
今回は以下のようなグラフを例にノードb（Start）からノードg（End）への探索を実装してみます。

入力

以下のようなDictionary形式で各ノードと隣接するノードを表すことにします。
search.py
graph = {
    'a': ['c'],
    'b': ['c','f'],
    'c': ['a','b','d'],
    'd': ['e','f'],
    'e': ['d'],
    'f': ['b','d','g'],
    'g': ['f']
}
深さ優先探索(Depth-First Search)

考え方

Start(b)の隣接ノードそれぞれについてEnd(g)と一致しているかチェック
【一致しない】→Startの隣接ノードをStartにして探索の関数を再帰する
【一致する】→終了

（mark_dictに既出のノードを記録して無限ループを回避している）

実装
search.py
import sys

# 深さ優先探索
def DFS(graph,start,end,mark_dict):
    mark_dict[start] = True
    for target in graph[start]:
        if target == end:
            print('Exist')
            sys.exit()
        elif target not in mark_dict:
            mark_dict[target] = True
            DFS(graph,target,end,mark_dict)

DFS(graph,'b','g',{})
幅優先探索(Breadth-First Search)

考え方

Start(b)の隣接ノードがEnd(g)と一致しているかチェック
【一致しない】→Queueにその隣接ノードに隣接するノードを入れ、Queueの最初の要素をStartにしてループさせる
【一致する】→終了

実装
search.py
# queueの実装
class Queue:
    def __init__(self):
        self.queue = []
    # 先頭にデータを追加
    def add(self,x):
        self.queue.append(x)
    # 最後のデータを削除
    def remove(self):
        del self.queue[0]
    # 最後のデータを返す
    def peek(self):
        return self.queue[0]
    # 空かどうか確認
    def is_empty(self):
        return len(self.queue) == 0

# 幅優先探索
def BFS(graph,start,end):
    mark_dict = {}
    queue = Queue()
    queue.add(start)
    while queue.is_empty() == False:
        target = queue.peek()
        queue.remove()
        if target not in mark_dict:
            mark_dict[target] = True
            if target != end:
                for elem in graph[target]:
                    queue.add(elem)
            else:
                print('Exist')
                sys.exit()

BFS(graph,'b','g')
双方向探索

幅優先探索を元のノードと目的地のノードの二ヶ所から同時に行う探索方法。
なぜこの探索方法が早いのかについての議論を書いておく。
各ノードがk個の隣接ノードを持っているグラフだと仮定する。
・幅優先探索だと、深さdの時点でk^dのノードを調べる必要がある。
・双方向探索だと、2つの探索がおよそ深さd/2の時点で衝突する。探索ノード数は元のノードからがk^(d/2)、最終ノードからk^(d/2)となり、合計ノード数は2*k^(d/2)となり幅優先探索よりも少なくて済む。

こちらの記事を参考にした。
https://qiita.com/guicho271828/items/cf1e7fcb98b1f074eacb

考え方

Start(b)由来かEnd(g)由来か識別するためのFlagを格納するQueueのClassを作る。
あとは幅優先探索とほぼ同じ考え方
Queueから要素取り出し
【自分の探索で未出】→Queueにその隣接ノードに隣接するノードを入れ、Queueの最初の要素をStartにしてループさせる
【自分の探索で既出】→Pass
【反対の探索で既出】→終了

（それぞれの探索で既出のノードを格納するmark_dictを設定した）

実装
search.py
# 双方向探索のためのフラグ付きqueueの実装
class FlagQueue:
    def __init__(self):
        self.queue = []
    # 先頭にデータを追加
    def add(self,x,flag):
        self.queue.append([x,flag])
    # 最後のデータを削除
    def remove(self):
        del self.queue[0]
    # 最後のデータを返す
    def peek(self):
        return self.queue[0]
    # queueをprint
    def get_queue(self):
        print(self.queue)
    # 空かどうか確認
    def is_empty(self):
        return len(self.queue) == 0

# 双方向探索
def BidirectionalSearch(graph,start,end):
    mark_dict = {0: {},1: {}}
    queue = FlagQueue()
    queue.add(start,0)
    queue.add(end,1)
    while queue.is_empty() == False:
        target = queue.peek()
        queue.remove()
        # 'target'が反対側のdictionaryで既出
        if target[0] in mark_dict[1-target[1]]:
            print('Exist')
            sys.exit()
        # 'target'が自分側のdictionaryで未出
        elif target[0] not in mark_dict[target[1]]:
            mark_dict[target[1]][target[0]] = True
            for elem in graph[target[0]]:
                queue.add(elem,target[1])


BidirectionalSearch(graph,'b','g')

　 Tweet
Qiitaで続きを読む

投稿日：2019-10-24T00:17:48+09:00

Slackにpythonからメッセージを送信する

概要

slackのapiをたたいてslackにpythonからメッセージを送信する。

やること

1. slackのワークスペースを作成

2. slack-appを作成

3. python-apiを実装

1. slackのワークスペースを作成

2. slack-appを作成

参考 (Slack API 推奨Tokenについて)

2-1. Slack apiにアクセスする。

2-2. install your app to your workspace

名前などを決めてappを作成する。

2-3. token取得

install Appの
「Bot User OAuth Access Token」を使用(botの発言)。
「OAuth Access Token」を使用すると個人の発言になる。

2-4. 権限の追加

OAuth PermissionsのScopesで以下の権限を追加する。
- chat:write:user (Send messages on the user’s behalf)

3. python-apiを実装
import requests


class SlackDriver:

    def __init__(self, _token):
        self._token = _token  # api_token
        self._headers = {'Content-Type': 'application/json'}

    def send_message(self, message, channel):
        params = {"token": self._token, "channel": channel, "text": message}

        r = requests.get('https://slack.com/api/chat.postMessage',
                          headers=self._headers,
                          params=params)
        print("return ", r.json())


if __name__ == '__main__':
    token = ''  # TODO your token.
    slack = SlackDriver(token)
    slack.send_message("Hello World! from python", "#random")
tokenをセットする。

参考

必要な Slack API はどれ？

Slack API 推奨Tokenについて

chat.postMessage api (機能を追加するならよく読む)

Slack api

https://github.com/st34-satoshi/slack-api

終わりに

python-slackclient というとても便利なライブラリがあるが、うまく使えなかった。sslがおかしいとか言われる...

これから

slackからのメッセージを受信できるようにしたい

　 Tweet
Qiitaで続きを読む

20191024のPythonに関する記事は30件です。

Jupyter Notebook で "NameError: name '__file__' is not defined" となる

疑問

原因

解決

参考

DeepChemでCross-Validationをやってみる

はじめに

環境

ソースコード

おわりに

参考

Zendeskのチケット検索で楽するためのツール simple-zendesk-searcher（仮称）を作った

結論

なぜ作ったのか

使い方 & ソース

まとめ

[ゼロから作るDeep Learning]勾配式実装準備 Pythonで微分と偏微分を実装してみた。

はじめに

微分の知識が必要な理由

微分

偏微分

事前分布に依存しないベイズ推定＆モデル比較

概要・導入

理論

例１：ガウス分布の混合比推定

例２：ガウス分布の混合比＆分散推定

サンプリング用コード

結果のプロット：

散布図：

Ralative belief updating function:

結局なんなのよ

まとめ

pylint をインストールしようとして sitecustomize のエラーに悩まされた話

概要

詳細

参考

Python よくわかるように書いた問題集no.1

問題1

出力結果

ソースコードの説明

変数の定義

フローチャート

大規模日本語ビジネスニュースコーパスを学習したXLNet（MeCab+Sentencepiece利用）モデルの紹介

はじめに

XLNetとは

弊社のXLNet事前学習済モデルの特徴

学習元データ

tokenizer

文字コード

事前学習済モデルの詳細

ダウンロード

利用方法

XLNetの日本語事前学習済モデルを利用した所感

BERTとの比較

NFKDとNFKCの比較

が→か 、で→て 等の助詞でのロスト

サービス→サーヒス 等の一般名詞や固有名詞でのロスト

まとめ

statsmodels:一般化線形モデル

一般化線形モデル(Generalized linear model)

例：教育の評価に関する分析(二項分布)

例：スコットランの選挙 (ガンマ分布)

例：非正準リンクをもつガウス分布

関連・参考サイト

参考

セミナー案内

flask_sqlalchemyで複数のデータベースに接続したい

ディレクトリ構成

コンフィグの書き方

モデルの書き方

printやdebugやpdbはなんかいや

参考

単純なテスト

最強のユニットテスト

TestLoaderの場合の出力

unittest.main()の場合の出力

特定のメソッドのみ実行の場合の出力

保有数が0の場合の売り注文を回避する #QuantX

order is ignored due to amount is 0

Jupyter Notebook で "NameError: name 'file' is not defined" となる

が→か、で→て等の助詞でのロスト

サービス→サーヒス等の一般名詞や固有名詞でのロスト

　単純なテスト

マルチスレッドは良い

Ⅰ.　ArcGISからQGISへ

Ⅱ-1.　QGISのPythonコンソールでの作業

Ⅱ-2.　ArcGISのPythonコンソールでの作業

8章　例外

9章　クラス　