StatsBeginner: 初学者の統計学習ノート

初学者が統計学、機械学習、R、Pythonの勉強の過程をメモっていくノート。

Python3でのメール送信時に日本語の差出人名を使う

以前のエントリで、Pythonからのメールの送信を試しましたが、


www.statsbeginner.net


この時は文中にも書いているとおり、差出人名を日本語表示するのがうまく出来ませんでした。
ところがその問題は、下記の方法で解決しました。


teratail.com


要するに、'差出人名 <アドレス>'という文字列を用意する時に、まるごとMIMEエンコードしてはダメで、「差出人名」のところだけがエンコードされるようにしなければならない。
それは前回のエントリ時点でも分かっていたんですが、なんか書き方が間違ってたようです。
どうやら、まず差出人名の文字列を.encode()メソッドでエンコードした上で、それをHeader()関数に与える際もHeader()関数の引数にエンコードを指定し、さらに、Headerオブジェクトそのものに.encode()メソッドを書き加えるということをしなければならなかったようです。
ただ、1個めの.encode()は要するに、Header()関数の引数としてエンコードされたバイト文字列を渡しているわけですが、べつにここにstr型で渡しても問題なく作動しました。
最後の.encode()は、これをつけることでRFCに沿ってMIMEエンコードされた文字列をstr型で受け取ることができるということのようです(説明)。これをつけないとこの部分がemail.header.Headerクラスのままになるので、あとで%sするときに不都合が起きるようです。


ということで、下記のように、

sender = '%s <%s>'%(Header('差出人名'.encode('iso-2022-jp'),'iso-2022-jp').encode(), メアド)


というように書いて、これをMIMEオブジェクトに与えて送信したら、ちゃんと表示されました。
一件落着です。


あと、何回か使ってみて、送信エラーが起きる理由は大きくわけて、

  1. 本文や名前(文中に宛名を差し込みする時)の日本語文字に'iso-2022-jp'でエンコードできないテキストが含まれている場合
  2. 送信途中でネットの回線が切れた場合

ですね。
なんかもうべつにUTF8で良い気はしてきました。

東京一極集中に関するデータ(ソースのメモ)

以下は、分析とかではなくデータのソースについてのメモである。
人口に関する「東京一極集中」がいかに凄いかは、世銀のサイトにある以下のデータをみるのが分かりやすい。


Population in the largest city (% of urban population) | Data


主要国では、日本の30%(都市人口の30%が東京圏に集中している)というのは極めて高いほうで、G20内ではアルゼンチンの35%負けるもののそれに次ぐ2位で、3位はサウジアラビアの20%。米・独・伊・中・露などは10%未満となっている。
ただ、後述するように韓国がソウルの人口しかカウントされておらず、首都圏という概念に広げると人口の半数を占めて堂々の1位となるようだ。


しかし韓国は国土面積が日本の4分の1ぐらいしかないし、サウジは砂漠だらけの国、アルゼンチンは歴史的にもブエノスアイレスという中世以来の貿易都市に後背地の草原をくっつけて国にしたようなところだから、日本のようにある程度広い国で歴史的にも有力な地方都市がある国の状況としては、異例ではある。しかも集中率が年々上昇している。(ドイツも年々上昇しているが、そもそもの水準が日本よりだいぶ低い。)


上記データはどういう由来かというと、国連の"World Urbanization Prospects"というレポートにおいて各国の「都市人口」が取りまとめられており、この数字を元に割合を出しているようだ。


World Urbanization Prospects - Population Division - United Nations


上記ページにまとまっているデータのうち、先ほどの世銀の集中率が使っているのは恐らく、

  1. 各国の、30万人以上が集積して住んでいる都市圏の、都市圏ごとの人口("Urban Agglomerations"中の"WUP2018-F12-Cities_Over_300K.xls")
  2. 各国の都市人口("Urban and Rural Populations"中の"WUP2018-F03-Urban_Population.xls")*1


の2つで、前者のうちその国で最大のものの人口を、後者で割っているのだと思う。分母が「都市人口」であって「総人口」ではないという点に注意が必要だ。日本の場合、その都市人口というのが1億1千万人ぐらいになっていたので、ほとんど全員が都市人口にカウントされているが。


なお、計算すると小数点以下が微妙に合わない。「年初」「年末」「年中」の違いとか、センサス未実施年のデータをどうするかとかの調整の都合が何かあるのかな?(具体的には確認していない。)
ちなみに、ここでいう「都市」というのは、"urban area"と"rural area"の対比でいう前者なので、「都会」と言ったほうが分かりやすいかもしれない。


都市(urban area)の定義は各国政府によるとされ、微妙な調整法などについては下記の2014年のレポートに注記がある。


https://esa.un.org/unpd/wup/publications/files/wup2014-highlights.pdf


日本の東京圏の定義がどうなっているかというと、先ほどの"WUP2018-F12-Cities_Over_300K.xls"のNOTES 70によれば、

(70) Major Metropolitan areas (M.M.A.) are defined by the Statistics Bureau of Japan. Census figures for 2005, 2010 and 2015 refer to the Kanto M.M.A.; figures from 1990 to 2000 are based on the Keihinyo M.M.A., and figures from 1960 to 1985 are based on the Keihin M.M.A. As a reference, the population of Tokyo-to was estimated at 12.1 million persons and of the Tokyo Ku-area at 8.1 million in 2000.


とされており、時期によって定義がちがうが、たとえば最新の"the Kanto M.M.A."とは、以下のページに載っている"Kantō Major Metropolitan Area (関東大都市圏)"であろう。


Greater Tokyo Area - Wikipedia


なお韓国のソウルについては「ソウル特別市」の人口が採用されている(NOTES 355)のだが、「首都圏」という意味ではもっと広いエリアが定義されるらしく、これだと総人口の半分ぐらいが首都圏に住んでいるらしいから、日本の東京集中度(約30%)を上回ることになるだろう。


首都圏 (韓国) - Wikipedia


また、台湾も、台北市が首都でその隣の新北市のほうが(恐らくベッドタウン化で)人口が多いのだが、この2つを合わせると集中度は38%ぐらいになる。ただし国土面積が日本の38万平方キロに対して、韓国は10万平方キロ、台湾はわずか3.6万平方キロと狭いので、日本とは事情が違うと考えたほうがよいだろう。


ところで、研究室に飾っておくために「立体日本地図」を買ったのだが、関東平野がいかに異常な広さかということがよく分かる。昔住んでいた茨城県つくば市では、目が良ければ筑波山の山頂から東京の新宿まで見通せるというのは有名な話だった。


f:id:midnightseminar:20180524173457j:plain


細かい定義を知らないが、関東平野の面積は1万7000平方キロメートルで、四国の1万8000平方キロメートルにほぼ匹敵する。
こんなに広い平地は北海道にすらなく、要するに日本の国土というのは、「殆どが山であり、点々と平野があってそこに都市が発達してきたが、関東にだけ四国と同じ広さの巨大平野がある」という構造になっている。


平野の定義がよくわからないのだが、以下のページに載っている都道府県別の地形・傾斜別面積をみると、関東1都6県の「丘陵地」「台地」「低地」面積を合わせると1万8,571平方キロメートルになる。


統計局ホームページ/第1章 国土・気象


日本全体では13万7,768平方キロだから、だいたい平野の13.5%が関東地方にあるということになる。
以下のページの「自然環境」のデータをみると都道府県別の可住地面積が分かるのだが、日本全体では12万4,038平方キロメートルであるところ、関東1都6県で1万8,261平方キロメートルなので、こちらは14.7%が関東にあるということになる。


社会・人口統計体系 統計でみる都道府県のすがた2018 | ファイルから探す | 統計データを探す | 政府統計の総合窓口


最初にあげた東京一極集中率の30%には、茨木・栃木・群馬などは含まれない*2はずだから、東京・千葉・神奈川・埼玉の1都3県の可住地面積が全国に占める割合を見てみると、だいたい7%ぐらいだ。
関東平野が恐ろしく広いとは言え、7%の可住地に30%の人が住んでいるのは、やはり偏っているようには思える。しかし関東平野全体ではまだ余力があるとも言えるので、北関東へのアクセスが劇的に改善した場合、さらに尋常ではない都市圏が形成されるという可能性もあるのかもしれない。

*1:30万人未満の「都市」人口も含まれている

*2:茨木南部は微妙に入る?未確認

RでVARモデル&インパルス応答関数を求めるとき、「ショックの値」をどうやって出すのか?

RでVARモデルを推定してインパルス応答関数を出す時に、インパルス応答関数が対応しているところの「ショック」の大きさが幾らなのかを、どこから得たらいいのかという疑問がわきました。


f:id:midnightseminar:20180504222113p:plain


こういう図を見せられたとして、まぁ増えるのか減るのかさえ分かればいい場合が多いかもしれませんが、タテ軸の数字はレスポンス側の変数の値だから意味が分かるものの、これがインパルス側の変数の当期の値に「どれだけのショック」を与えた場合に生まれる変動なのか数字で言ってくれ、って頼まれたらどう答えるのかなと。


で、{vars}パッケージの仕様書を見ると、コレスキー分解をしているようなので、たぶんショック側の変数の「撹乱項の標準偏差」だろうと思うのですが、確かめ方が私には俄かには分からなかったので、不躾ながらパッケージの開発者(Bernhard Pfaff氏)にメールで「ショックの大きさは、撹乱項の標準偏差の推定値ってことで合ってますか?」と聞いてみました。
そしたら親切にも速攻で返事が返ってきて、「そのとおりや! "vars:::Psi.varest" と打って俺のPsiのコードを見とけ!」と言われました。
ここでいうPsiは、


https://cran.r-project.org/web/packages/vars/vars.pdf


この仕様書の28ページの上の方に載ってる式のΨのことで、この行列の中身が直交化インパルス応答に対応しています。なおこれはVARのVMA表現ですね。*1
で、コードをみてみると、

> vars:::Psi.varest
function (x, nstep = 10, ...) 
{
    if (!(class(x) == "varest")) {
        stop("\nPlease provide an object of class 'varest', generated by 'VAR()'.\n")
    }
    nstep <- abs(as.integer(nstep))
    Phi <- Phi(x, nstep = nstep)
    Psi <- array(0, dim = dim(Phi))
    params <- ncol(x$datamat[, -c(1:x$K)])
    sigma.u <- crossprod(resid(x))/(x$obs - params)
    P <- t(chol(sigma.u))
    dim3 <- dim(Phi)[3]
    for (i in 1:dim3) {
        Psi[, , i] <- Phi[, , i] %*% P
    }
    return(Psi)
}


ここでsigma.uとなっているものが、残差の分散共分散行列ですね!
ということは、この行列と同じものを自分で求めて、対角要素の平方根を取ればショックの値になるはずです。
これは、VARモデルの推定結果であるx(varestクラスオブジェクト)の中に入っている残差の行列


resid(x)


を取り出してクロス積をとって自由度で割ったものです。
自由度は、元のデータがN変量×T期分あって、VAR(p)モデルを構築したのであれば、


df = (T-p) - (Np + 1) # +1は定数項の分


になります。
以下、念のため実際にサンプルデータで推定してみてやり方を確認しておくことにします。


まずデータを読み込みます。変数が4個ありますが、労働生産性と実質賃金だけ使おうと思います。

> library(vars)
> 
> # 練習データのCanadaを時系列データ型でインポート
> # e: 千人あたりの雇用
> # pros: 労働生産性
> # rw: 実質賃金
> # U: 失業率
> 
> data(Canada)
> 
> # 私はts型に慣れてないのでデータフレームにしますw
> canada <- as.data.frame(as.matrix(Canada))
> 


どんなデータか描画しておきます。

> split.screen(c(2,1))  # 描画デバイス分割
> screen(1)
> plot(canada$prod, type='l', main='Productivity')  # 労働生産性
> screen(2)
> plot(canada$rw, type='l', main='Real Wage')  # 実質賃金


f:id:midnightseminar:20180504222758p:plain


では、VARモデルを推定していきます。
簡単にするため、労働生産性(prod)と実質賃金(rw)という2つの変数だけでやってみます。

> ### 労働生産性と実質賃金のデータだけでVARモデルを構築します
> # ラグの選択
> lag <- VARselect(canada[,c(2:3)], lag.max=10)$selection[1]  # AIC最適ラグ数を選択
> print(lag)
AIC(n) 
     3 
> 
> # VARの推定
> var.canada <- VAR(canada[,c(2:3)], p=lag, type='const')
> summary(var.canada)

VAR Estimation Results:
========================= 
Endogenous variables: prod, rw 
Deterministic variables: const 
Sample size: 81 
Log Likelihood: -174.944 
Roots of the characteristic polynomial:
0.9843 0.813 0.813 0.5431 0.5431 0.4038
Call:
VAR(y = canada[, c(2:3)], p = lag, type = "const")


Estimation results for equation prod: 
===================================== 
prod = prod.l1 + rw.l1 + prod.l2 + rw.l2 + prod.l3 + rw.l3 + const 

        Estimate Std. Error t value Pr(>|t|)    
prod.l1  1.21318    0.11423  10.620   <2e-16 ***
rw.l1   -0.02120    0.09222  -0.230   0.8188    
prod.l2 -0.21853    0.18092  -1.208   0.2309    
rw.l2   -0.14257    0.13835  -1.031   0.3061    
prod.l3 -0.04157    0.11778  -0.353   0.7251    
rw.l3    0.16883    0.08846   1.909   0.0602 .  
const   17.21994   10.69523   1.610   0.1116    
---
Signif. codes:  
0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 0.6841 on 74 degrees of freedom
Multiple R-Squared: 0.976,	Adjusted R-squared: 0.974 
F-statistic: 501.3 on 6 and 74 DF,  p-value: < 2.2e-16 


Estimation results for equation rw: 
=================================== 
rw = prod.l1 + rw.l1 + prod.l2 + rw.l2 + prod.l3 + rw.l3 + const 

        Estimate Std. Error t value Pr(>|t|)    
prod.l1 -0.12808    0.13563  -0.944  0.34808    
rw.l1    1.08616    0.10949   9.921  3.1e-15 ***
prod.l2 -0.27601    0.21480  -1.285  0.20281    
rw.l2   -0.17954    0.16426  -1.093  0.27793    
prod.l3  0.44228    0.13984   3.163  0.00227 ** 
rw.l3    0.06736    0.10502   0.641  0.52323    
const   -3.05700   12.69828  -0.241  0.81042    
---
Signif. codes:  
0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 0.8122 on 74 degrees of freedom
Multiple R-Squared: 0.9987,	Adjusted R-squared: 0.9985 
F-statistic:  9176 on 6 and 74 DF,  p-value: < 2.2e-16 



Covariance matrix of residuals:
         prod       rw
prod 0.467974 0.002593
rw   0.002593 0.659676

Correlation matrix of residuals:
         prod       rw
prod 1.000000 0.004666
rw   0.004666 1.000000

> 


VAR(3)の推定が終わりました。
次に「労働生産性→実質賃金」という方向の、直交インパルス応答関数を求めます。あわせてグラフの描画もしておきます。

> irf.canada <- irf(var.canada,impulse='prod', 
+                   response='rw', 
+                   ortho = TRUE, 
+                   n.ahead=10,ci=0.95,
+                   cumulative = FALSE,
+                   runs=300
+                   )
> 
> print(irf.canada$irf)
$prod
                rw
 [1,]  0.003789932
 [2,] -0.083499262
 [3,] -0.386475239
 [4,] -0.440968850
 [5,] -0.394311432
 [6,] -0.341204111
 [7,] -0.296874634
 [8,] -0.234418701
 [9,] -0.160965629
[10,] -0.089254559
[11,] -0.023464820

> plot(irf.canada)


f:id:midnightseminar:20180504222113p:plain


これでつまり、当期(第1期)のprod(労働生産性)に撹乱項の1標準偏差分のショックが加えられた時の、各期のrw(実質賃金)の増減が得られました。


で、問題は「その労働生産性の撹乱項1標準偏差分ってどんだけやねん」ということですが、これは次のようにして計算できます。回りくどい書き方ですが、

> # 自由度
> Obs <- var.canada$obs  # 推定されるyの期数(元データの行数-ラグ次数)
> N  <- var.canada$K     # 変数の数
> p <- lag               # ラグ次数
> df <- Obs - (N*p + 1)
> 
> # 残差の分散共分散行列
> sigma.u <- crossprod(resid(var.canada))/df
> print(sigma.u)
            prod          rw
prod 0.467973610 0.002592639
rw   0.002592639 0.659676363
> 
> # irf()関数でそれぞれをimpulsに設定したときのショック=撹乱項の標準偏差
> shock.prod <- sqrt(sigma.u[1,1])
> shock.rw <- sqrt(sigma.u[2,2])
> 
> print(shock.prod)  # prodの撹乱項1標準偏差分の値
[1] 0.684086
> 


ここまでで、ショックの値は計算できました。
つまりさっきのインパルス応答関数と合わせていうと、当期の労働生産性の撹乱項に0.684086のショックを加えると、当期の実質賃金には0.003789932、2期目には-0.083499262...の影響があるということが分かったわけです。(実際は先に単位根とか共和分とかを見ないといけませんがここでは分析自体が目的ではないので省略。)


ちなみにこの撹乱項の標準偏差の値、どっかに入ってるのではないかと思って探してみたら、いちいち計算しなくても、以下のようにすれば一発で取り出すことができました。実務上はこれを使えばいいと思います*2。summaryをprintしたときの出力にも、'residual standard error'として載っています*3
念のためですが、さっきのsigma.uの対角要素は分散、こっちのsigmaは標準偏差の推定値です。

> summary(var.canada)$varresult$prod$sigma
[1] 0.684086
> 
> summary(var.canada)$varresult$rw$sigma
[1] 0.8122046


念のため、これがインパルス応答関数のショックの値に使われているという理解であってるのかどうか、


沖本(2010).経済・ファイナンスデータの計量時系列分析,朝倉書店.


に書いてある方法に照らしながら手計算でインパルス応答を出してみて、それが{vars}パッケージのirf()関数の結果と確認するかを見ておきます。ただし後述の通り、使われているショックの値が合ってるかどうかを確認するだけなら、1年後、2年後…の逐次計算をしなくても、当期のresponseだけ見れば分かります。


まず、沖本(2010)のp.87-90あたりに書いてある「三角分解」(ほかの呼び方もあるようですが)の行列Aは、以下のように作れますね。

a <- sigma.u[2,1]/sigma.u[1,1]
A <- matrix(c(1,a,0,1),2,2)


で、1つ目の変数から2つ目の変数*4へのインパルス応答のうち、「当期の値」(今回の例でいうと1期目のrwのレスポンスの値)については、「2つ目の変数の撹乱項のうち、1つ目の変数から影響を受けてる部分」がその値に一致すると思います。
沖本(2010)のp.89の、


f:id:midnightseminar:20180504223144p:plain


この部分に「0.3u_1t」ってのがありますが、要するに2つ目の変数の撹乱項ってのは、2つ目の変数に固有の撹乱項と、1つ目の変数の撹乱項を0.3倍した成分に分かれるということですね。これを分けたのが直行化ってやつです。
で、インパルス応答関数とVARモデルの定義から言って、2つ目の変数の1期目のレスポンスには、1つ目の変数の1期目に与えられたショックを0.3倍(今回の場合でいうとa倍)した値だけが入ることになるので、これだけなら手計算がすぐできます。2期目以降の手計算は面倒ですが。
今回の例でいうと、以下のようになります。

> response.rw.1 <- a * shock.prod
> print(response.rw.1)
[1] 0.003789932
> 


この通り、さっきの「0.003789932」と一致してるので、合ってると思います。
たぶん・・・。
(時系列データの分析は実務でほとんどやらないので、何か根本的に間違ってたらごめんなさい。)

*1:VMA表現についての解説はハミルトン『時系列解析』の10-11章を見るのがいいと思います。

*2:個人的には、撹乱項の標準偏差を推定する上で自由度をいくらにするか自信がなかったので、確認できてよかったですw

*3:eとかrwの推定値の標準誤差だから、撹乱項の標準偏差のことになります。

*4:沖本(2010)のp.90でいう再帰性、つまり変数の順序性があるので、1つ目、2つ目という意識は意外と重要ですね。

Rでの単位根検定はadf.test()関数よりCADFtest()関数がいいのでは?

時系列データをあまり扱わないのでまじめに考えてなかったんですが、Rで単位根検定をする場合、拡張ディッキー=フラー検定(augmented Dickey–Fuller test)を実施してくれるadf.test()という関数があります。
しかしこの関数は、

  • 考慮するラグの次数を指定しなかった場合、サンプルサイズ(時系列データの長さ)を基準にして自動選択している
  • 定数項もトレンド項も含むパターンしかやってくれない


という制限があり、後者についてはたいていの場合含めといたほうがいいらしいので問題ないとして、前者はよくわからない基準であって*1、本来はAICなどに基づいて選択されたラグ次数を用いたモデルで検定したほうがいいと思われます。実際ラグ次数によって結果がけっこう変わることもあるので、慎重に選んだほうがいいかと思います。


ということで、以下の記事で説明されているように、CADFtestパッケージのCADFtest()関数をつかっとけば簡単でいいんじゃないでしょうか。


定常過程かどうかのチェック(ADF検定)


上の記事だと色々やってて分かりづらいかもしれないので、パッケージの仕様書を見たほうがいいと思いますが、


Package ‘CADFtest’


これを読むと、共変量付きのADF検定ってのもできるらしいです。てうかそもそも,頭の"C"がCovariateのCで、それ用のパッケージみたいですね。
で、一番単純な使い方としては、

#dに時系列データがベクトルで入ってるとする

CADFtest(d, 
   type="trend",   # トレンド項も定数項もあり
   max.lag.y=10,  # ラグの最大次数を自分で適当に指定
   criterion='AIC'  # ラグ次数はAIC基準で選ぶ
)


でいいんじゃないでしょうか(共変量を考えるときは、dを入れている最初の引数にモデル式を記述するとのこと)。
typeのところを"trend"にすると定数項もトレンド項もあり、"drift"にするとトレンド項なし、"none"にするとどっちもなしになるようです。よほどの理論的理由がない限り、両方入れておくのが保守的らしいです。
単純に結果だけ知りたい場合は、上記をprintして出てくるp-valueが0.05を下回っていれば、5%水準で「単位根あり」の帰無仮説を棄却できることになります。0.05より大きければ、単位根過程の疑いありです。


詳しい使用例が以下のドキュメントにのってました。
Covariate Augmented Dickey-Fuller Tests with R


間違ってたらごめんなさい。

*1:「このRの自動選択基準はきわめていい加減なので、信じてはいけません」と評している方もいます(リンク

Rで要素番号の指定の仕方をミスった

 考えてみればそりゃそうか、という感じではあるのですが、またいつかミスりそうなのでメモしておきます。
 たとえば以下のような感じで、startとendの値を変えて適切な期間を取りたいとします。

> v <- c(1961, 1962, 1963, 1964, 1965, 1966, 1967, 1968, 1969, 1970)
> start <- 3
> end <- 5
> v[start:end]
[1] 1963 1964 1965


 それで何かの都合から、startとendの値を少しずらしたいことがあったとします。
 たとえば終了年度を1年早くして、1963年から1964年を取りたくなって、end-1みたいなことをするとミスります

> v[start:end-1]
[1] 1962 1963 1964


 このように、意図としては(1963, 1964)という出力が欲しかったのに、(1962, 1963, 1964)になってしまっています。
 なんでこうなるかというと、上の書き方だと、start:endの部分でまず(3, 4, 5)というベクトルが生成され、この各要素から1を引くという処理になってしまって(2, 3, 4)になるからですね。
 ちなみに対処としては、start:c(end-1)にしておけば意図どおりになります

 同様の理屈で、

> v[start+1:end]
[1] 1964 1965 1966 1967 1968


 これは1:endつまり(1, 2, 3, 4, 5)のそれぞれの要素に、start=3を加えるという処理になっています。

> v[start-1:end]
Error in v[start - 1:end] : 
  only 0's may be mixed with negative subscripts


 これはエラーになっていますが、なぜかというと、1:endの部分で(1, 2, 3, 4, 5)という数列が生成され、start=3からそれぞれを引いた数の列になるからです。
 すなわちこの場合、(2, 1, 0, -1, -2)という数列になり、要素番号が適切に指定できなくなります。

Pythonの簡単なコードでメールを自動送信してみる

意外と簡単にできた

 メールを300人ぐらいに発信する必要がありまして、Toに全員入れるわけにはいかないし、BCCで送るのもダサいかなと思って、「1人1人を個別にToに指定して、同じ件名・同じ文面のメールを送る」ってのをPythonでやってみました。*1
 標準モジュールのemailってのとsmtplibってのを使って、50行程度のコードで簡単に送れました。
 1点心残りなのは、後述のとおりFromの欄に日本語の差出人名を表示させるやつが、色々調べたものの結局できませんでした。

用意するもの

 アドレスリストをCSVで用意して、本文はテキストファイルに書いておきました。

f:id:midnightseminar:20180212202431p:plain
f:id:midnightseminar:20180212202409p:plain

コード

 ネットでsmtplibを使ったPythonでのメール送信の解説を探すと、sendmail()というメソッドで送っているものと、send_message()というメソッドで送っているものがあります。sendmail()が基本なのですが、send_message()はそれをより簡単に扱えるようにラップしてくれているような感じらしいです。今回はsend_message()を使いました。


 コードは以下の通りですが、全体としては、

  1. モジュールを読み込む
  2. 差出人アドレスやメールサーバの認証情報などの設定項目をまとめて書いておく
  3. メールサーバに接続する
  4. emailモジュールのMIMETextでテキスト形式のメール本体(MIME文書)をつくる
  5. MIME文書に、差出人、件名、宛先等の情報を入れていく
  6. smtplibのsend_messageでメールを送信する
  7. サーバとの接続を終了する

 という流れになっています。

# モジュールの読み込み
import time
import smtplib
from email.mime.text import MIMEText
from email.header import Header
import pandas as pd

# 基本的な設定たち
srv_smtp = 'XXXXXX.XXXXXX.jp'  # SMTPサーバ
srv_port = 587                 # ポート番号
srv_user = 'XXXXXX'            # サーバのユーザ名(ちなみに私が使ってるやつだとメアドがユーザ名)
srv_pw   = 'XXXXXXXX'          # サーバのパスワード
jp_encoding = 'iso-2022-jp'    # 日本語文字エンコーディングの指定
add_sender = 'XXXX@XXXXXX.jp'  # 差出人(自分)アドレスの設定
add_bcc = 'XXXX@XXXXXX.jp'     # BCCの複製を送るアドレス
add_rcp_path = '/XXXXX/XXXXX/address_test.csv'  #アドレス一覧が入ったCSVの置き場
body_path =    '/XXXXX/XXXXX/body_test.txt'     # 本文を書いたテキストファイルの置き場
mail_subject = 'くさめの件につきまして'         # 共通の件名

# 本文ファイルの読み込み
with open(body_path, 'r', encoding='utf-8') as file:
    mail_body = file.read()


# 宛先リスト読み込み
# 元ファイルには名前も入れてるけどとりあえず使わないことにする
add_rcp_df = pd.read_csv(add_rcp_path, encoding='utf-8')  # csvの読み込みはpandasでしかやったことないので…
add_rcp_list = add_rcp_df['Address'].tolist()


# SMTPサーバへの接続
server = smtplib.SMTP(srv_smtp, srv_port)
server.ehlo()
server.starttls()  # TLSでアクセス
server.ehlo()
server.login(srv_user,srv_pw)  # ログイン認証

# 送信を繰り返す
for add in add_rcp_list:
    try:
        msg = MIMEText(mail_body.encode(jp_encoding), 'plain', jp_encoding,)
        msg['From'] = add_sender
        msg['Subject'] = Header(mail_subject, jp_encoding)
        msg['Bcc'] = add_bcc
        msg['To'] = add
        server.send_message(msg)  # 送信する
        time.sleep(3)  # 3秒まつ
    except:
        # なんかあった時用
        print('An error occured when sending a mail to ' + add)
        

# サーバ接続を終了
server.close()


 MIMETextでMIME文書を作って、件名等の情報を入れた後、特定の情報だけ差し替えるというのが上手く行かなかったので、forループでメールを1通1通送る際に、MIME文書の生成自体をまるごとやり直しています。これが適切なのかよく分かってませんがとりあえずメール送信には成功しました。
 あと、メールソフトに送信済みメールが残らないので、記録用にBCCで自分のアドレス宛にメールを飛ばしています。
 1万件とか送るのであれば、受信側のメールサーバにSPAM判定されないように時間を空けて送る必要があるかと思いますが*2、300件ぐらいならべつに全部即時送信してしまっていいような気はします。上記では一応3秒ずつ空けています。

なんか止まってた

 300件の送信中、2回止まりました。
 たぶん回線が不安な環境だったので、ネット接続が切れたんだと思いますw
 エラーが出たアドレスのところを確認して、あと一応BCCでコピーを飛ばしていたメールも確認して、未送信の人だけのアドレスリストを作ってやり直しました。
 上記のコードでは、何らかのエラーが出た時の対処はまったく記述していないので、簡単なコードでむやみに大量のメール送信をすると何が起きるかわからんという点には注意が必要かと思います。

日本語文字のエンコード

 エンコーディングのところに関して、参考にしたブログ記事などが一様に、日本語のメールで伝統的に使われているという'iso-2022-jp'を指定していたので、上記コードではその通りにしてますが、↓のページに書かれているように、今はべつにUTF-8でも問題ないようです。
 実際、UTF-8でも自分の持っている幾つかのメールボックスに送ってみましたが、ちゃんと見れました(ただ、iPhoneとMacでしか確認してないので見れて当然なのかもしれません)。
 
Pythonで日本語メールを送る方法をいろいろ試した
 
 受信側の環境にも拠るのかもしれないので、私は念のため伝統的なほうを使いましたが、実際どっちのほうが安全なのかはよく知りません。
 
 

差出人表示を日本語でする方法が分からない

 上記コードで、

add_sender = 'XXXX@XXXXXX.jp'

のところを

add_sender = 'K.Yoshida <XXXX@XXXXXX.jp>'

にすると、受信側のメールソフトで差出人名を「K.Yoshida」として表示してくれたりします。


f:id:midnightseminar:20180212202526p:plain:w300


 それは今回もできたんですが、ここに日本語の名前を入れる方法というのが色々難しく、結局ちゃんとはできませんでした。
 下記のようなページで紹介されているように、

  • send_message()ではなくsendmail()で送る
  • アドレス部分はそのままに、日本語の差出人名の部分だけエンコードする。その際単に文字列としてではなくHeaderインスタンスとして生成する

 という方法でできるらしいのですが、自分でやってみたところ、受信側サーバによっては受信拒否、一応送れたサーバでも、差出人名のあとに「@」と受信サーバの情報を表す文字列がくっついた状態で表示されてしまい、要するに不正なメールとして判定されたんだと思います。この辺、正式にできるやり方を調べる必要があります。
 Stackoverflow等をみると海外でもウムラウト付きの文字などを表示させようとして苦労している人がいました。


Python3 日本語でメール送信 - textbook
Pythonでメール送信時に送信者に日本語を使用する : fujishinko 雑記帳


[追記]
下記のとおり、成功しました。
Python3でのメール送信時に日本語の差出人名を使う - StatsBeginner: 初学者の統計学習ノート
[/追記]

 

拡張

 上記のようにとりあえず単純なメールを送れる状態にしておけば、あとは

  • アドレスリストから名前を取得して、本文の一行目に「◯◯様」と可変で入れる。
  • そもそも本文自体を人に合わせて変える。
  • 途中でネット回線が切れたりするのが怖いのでAWS等の仮想マシンから発射する。
  • multipartして、base64した添付ファイルをつける


 など、今後いろいろ工夫できるかと思いました。
 ただ、私のような素人が自前のコードでメール送信をすると、なんか事故が起きそう(宛先と添付ファイルが1個ずつズレるとか)なのでなるべくやりたくはないですね。
 
 

関連リンク

mimeTEXTの説明書
https://docs.python.jp/3.3/library/email.mime.html

smtplibの説明書
https://docs.python.jp/3/library/smtplib.html

ここにemailモジュールとsmtplibモジュールを使ってメールを送る際のコード例が載っています。
https://docs.python.jp/3.3/library/email-examples.html

ここに書かれているように、From欄とかnon-ascii文字を使いたければ、Headerモジュールを使ってHeaderインスタンスとして投入する必要があると書いてあります。
https://docs.python.jp/3.3/library/email.header.html

*1:MacのAutomatorでも、group mailっていうアクションを使うとアドレス帳で指定したグループあてに「1人1人をToに指定した別々のメール」として一斉送信ができるけど、冒頭に'Dear Bob,'みたいなgreetingを入れる必要があり、しかもこれが英語仕様しかないので強制的に最後にカンマが入るという、日本人にとっては中途半端な仕様です。

*2:1万件あれば、docomo.ne.jp等に同時に数百から数千通発射することになる

Macでの年賀状作成環境について(2018年版)

 ブログの趣旨と全然違いますが、備忘のためにまとめておこうと思います。
 私は、サボる年もありますが、出す時は年賀状を200人以上に出すので、けっこうな大仕事になっています。
 しかも東京周辺の知り合いが多く、3〜4年たつと半分ぐらいの人が引っ越してしまっている(自分と同年齢から下10年ぐらいまでの知り合いは、結婚によって引っ越すケースも多いので)イメージで、住所録のメンテナンスがけっこう大変です。

宛名職人にサヨウナラ

 昔は、年賀状は「宛名職人」(リンク)というMac用の年賀状ソフトで作成していました。宛名職人の住所録は、iCloudのアドレス帳に登録してある情報をvCardで書き出しておけば、差分を検出してどっちかで上書きするってのができて便利だったんですよね。
 私は一時期、年賀状ソフト、携帯電話、パソコンのメーラーなどに電話帳やアドレス帳が分散しているのが気持ち悪かったので全部統合してiCloudにぶち込んだことがあります。これで一元化されたので、iPhoneでアドレス帳を開くと、年賀状のやり取りがある人については住所も入っているという状態になりました。で、

  • 年末年始以外の時期に、誰かが引っ越したことを知った時は、iPhoneかMacでiCloud住所録を上書きする。
  • 年末に年賀状を準備する前に、iCloud住所録→宛名職人住所録の方向で上書きする。
  • 年賀状を出す。
  • 何人かは「あて所に尋ねあたりません」で戻ってきたり、向こうからくる年賀状によって私の認識している住所が異なることが判明する。
  • 届いた年賀状を見ながら、宛名職人の住所録画面を見ながら、カタカタ打ち込んで更新する。
  • 年賀状シーズンが終わったら、宛名職人住所録→iCloud住所録の方向で上書きする。


 みたいな感じでやっていました。
 ところが宛名職人は、毎年新しいバージョンのソフトを数千円で買わされるんですよね。古いバージョンでも動くときがありますが、Mac OSの新バージョンが毎年のように出るので、古いバージョンだと不具合が出るようになりました。
 しかも、何年前だったか忘れましたが、最新バージョンにしたのにトラブルが起きたことがあり、まぁ自分に原因があったのかも知れませんし詳細は忘れましたが、とにかく自分が安定して使いこなせないようなソフトは使うのをやめようと思って決別しました。

はがきデザインキットは使いこなせなかった

 その後2年ぐらいは、Windows機の筆まめを使ってやっていましたが、今年は「はがきデザインキット」(リンク)を使ってみようと思いました。Mac版があるので。
 郵政関係者に聞いたところ、はがきデザインキットは日本国内で最もたくさん利用されている年賀状ソフトらしいです。たしかに無料でこれだけの機能があれば、有料の筆まめとかを使う必要はあまりない気がします。


 ところが実際にトライしてみると、住所録の読み込みで躓きました。はがきデザインキット専用のcsvフォーマットなるものが配布されているので、それに従ってcsvを作成したのですが、文字コードを色々変えてみても読み込み後に文字化けしたので、4、5回トライして諦めました。


f:id:midnightseminar:20180103165708p:plain

宛名印字はプリントマジックで

 で、結論として今回は、宛名の印字にプリントマジック(リンク)というのを使うと、比較的簡単に作成できました。


f:id:midnightseminar:20180103165458p:plain


 こんな感じの画面で、CSVを読み込む際に、CSVのどのカラムをプリントマジックのどの項目に割り当てるかを選択します。
 ちなみに、元のCSVにはヘッダ行があるのでそれが1行目に表示されてしまっています。「一行目を無視する」をチェックすると消えますが、ヒモ付け作業が終わるまでは出しておいたほうがわかりやすいです。「項目3は姓に紐付ければいいんだな」とか分かるので。


 元の住所データを、どこで区切って保存しておくかというのは、宛名印字界隈では悩ましい問題ですね。そもそも日本の住所の決め方自体もけっこうややこしいので(以前エントリとしてまとめました)、どこで区切るのが正しいかといった定説もない気がします。


 私はもともと、


(1) 都道府県
(2) 市区町村(政令指定都市の場合は行政区まで入れる)
(3) 町・字から住居番号まで
(4) 建物名と部屋番号


 というふうに分けたデータを作っていましたが、「◯◯町1−2−3−405」みたいに、マンション名を省略してかつ部屋番号(ここでは405)をそのまま続けている住所しか把握していない場合は、(3)に「405」まで含めてしまったりしています。
 そもそも数字のところは「丁目」「番地」「号」という3点セットになっているとは限らないので、住居番号と部屋番号の境目は判別がしづらいです。


 で、今回は、プリントマジックでの印刷レイアウトを考えると、

  • 都道府県&市区町村
  • 町・字から住居番号まで
  • 建物名と部屋番号

 という3項目に分けたCSVを読み込ませると、いい感じでした。
 プリントマジックでは、たとえば項目9〜項目11をすべて「自宅住所」として指定すると、住所録画面ではこの3項目を連結した情報が表示されますが、印刷時にはこの3つを行で分けて印刷してくれます。
 その際、一行の長さをいい感じにするには、上記のような分け方がちょうど良かったです。


 住所録のうち、選択した行だけ印刷することもできるし、印刷するかしないかというデータ項目もあります。
 印刷しない人のデータなんてそもそも必要なの?と思われるかもしれませんが、「普段は出しているけど、今年はこの人から喪中のお知らせが届いたので出さない」みたいなパターンがありますし、こっちから出してないけど向こうからは届いたので年明けに「返事」を出すという場合は、その人たちだけ「印刷する」に設定して印刷するのが良いですね。


 プリンタの設定でハガキを指定すれば、ふつうに綺麗に印刷できました。印刷ソフトによっては、郵便番号の位置がプリンタとの相性で微妙にずれるため、微調整機能がついてたりしますが、今回は相性が良かったのか、何の問題もなく印刷できました。

デザイン面はPower Point

 宛名職人でも筆まめでもはがきデザインキットでもそうですが、テキストや画像を配置する操作に制限が多くてけっこうイライラしますよね。しかも年1回か(暑中見舞いを入れて)2回しか使わないので、UIのクセに習熟することがない。
 「だから全部フォトショで作る」という知り合いもいます。
 しかしフォトショなどの画像処理ソフトも、文章を多めに書く年賀状を作る場合は、テキストをいれる作業でちょっとイライラすると思います。また、文字も全部画像データとして吐き出されるというのはなんか気持ち悪い気もします。クッキリ印字するためには無駄に高解像度にしないといけないのではないか、とか。


 そこで今回は、Power Point(for Mac)で作ることにしましたが、これは大正解でした。写真や図形や文字の配置が自由自在で、とにかく使い慣れたUIなのでイライラが少ない!

  • ページ設定で、画面をはがきのサイズにする。
  • 写真や文字を配置する。この際、余白を考えずにパワポのページの端まで使う。
  • PDFとして保存する(余白のないPDFができる)
  • プリンタ(私はエプソンのもの)の設定で、余白3mmのハガキ印刷を指定。


 これでとてもいい感じに印刷できました。PDFなのでテキスト部分を無駄に解像度の高い画像にする必要もなし。


 来年もこのパターンでいくと思います。
 iCloudの住所録と同期するのは、もう色々めんどうなのでやめることにします。年賀状用のCSVを、それはそれとしてメンテナンスしていく所存です。