pandasでクロス集計

テーブルデータを分析する時、2つの列の値の関係を調べたいことが良くあります。
対象の列がどちらも連続値をとるのであれば、散布図を書きますが、
これが連続値ではなく性別や都道府県、世代、アンケートのn段階評価などの時は、散布図はあまり適しません。

このような場合は、クロス集計を使うと便利です。
(wikipediaでは分割表と呼んでいるようです。)

pandasにクロス集計を行う専用の関数があるのでそれを紹介します。
(職場ではcsvに書き出してtableauに読み込ませてやることも多いのですが、pythonで完結できればそれはそれ便利です。)

ドキュメントはこちら。
pandas.crosstab

aggfuncに値を渡せば、数え上げ以外にも色々できることがあり、pivot_table的な使い方もできるのですが、
とりあえずダミーデータを用意してデータの数え上げでクロス表を作ってみましょう。
ABテストなどで頻繁に使いますからね。


import pandas as pd
import numpy as np
# ダミーデータ生成
df = pd.DataFrame(
        {
            "data1": np.random.randint(1, 6, 100),
            "data2": np.random.choice(["foo", "bar"], 100),
        }
    )
# サンプル(5件)
print(df.sample(5))
# 出力
'''
    data1 data2
71      2   bar
13      1   foo
7       5   bar
16      1   bar
56      2   foo
'''

# クロス集計
cross_df = pd.crosstab(df.data2, df.data1)
print(cross_df)
# 出力
'''
data1   1   2   3   4  5
data2                   
bar    10  13   8  12  8
foo    12   9  11   9  8
'''

とても楽にクロス集計ができました。
SQLでこれを実装しようとするとかなりの手間なので、地味に重宝しています。

完全に蛇足ですが、crosstabを知らなかった頃はこういうコードを書いていました。
indexとcolumnsだけ指定して値が全部0のデータフレームを作って、
元のデータの値を順番に数え上げています。
今思えばcrosstab知らなくても、もう少しまともな書き方がありそうです。


cross_df_2 = pd.DataFrame(
        0,
        index=set(df.data2),
        columns=set(df.data1),
    )
for _, row in df.iterrows():
    cross_df_2.loc[row.data2, row.data1] += 1    
print(cross_df_2)
# 出力
'''
     1   2   3   4  5
foo  12   9  11   9  8
bar  10  13   8  12  8
'''

pythonで一般化線形モデル(GLM) – ポアソン回帰

以前の記事で、statsmodelsを使って線形回帰をやったので、今回はポアソン回帰をやってみます。
参考:statsmodelsで重回帰分析

データは久保拓弥先生の、データ解析のための統計モデリング入門 (通称緑本)の第3章から拝借し、
本に載っているのと同じ結果を得ることを目指します。
ちなみに本のコードはRで実装されています。
Rは勉強中なので写経はしましたがそれはそれとして、
今回はいつも使っているpythonのstatsmodelsでやってみます。

データはこちらのサポートページからダウンロードできます。
生態学のデータ解析 – 本/データ解析のための統計モデリング入門
3章の data3a.csv を保存しておきましょう。

このデータは、架空の植物の種子の数に関するもので、
種子の数が$y_i$, 植物のサイズが$x_i$, 施肥処理を行ったかどうかが$f_i$列に格納されています。

そして、種子の個数$y_i$が期待値 $\lambda_i=\exp(\beta_1+\beta_2x_i+\beta_3f)$のポアソン分布に従うと仮定した場合の尤度を最大化するように係数を決定します。
($f$は施肥処理を行ったら1, 行ってない場合は0)

早速やってみましょう。
ドキュメントはこのあたりが参考になります。
statsmodels.genmod.generalized_linear_model.GLM
Model Family に Poisson を指定すると、リンク関数は自動的にlogが選ばれます。


import pandas as pd
import statsmodels.api as sm

# データの読み込み
# ファイルの取得元
# http://hosho.ees.hokudai.ac.jp/~kubo/ce/IwanamiBook.html
df = pd.read_csv("./data3a.csv")

# f 列の値を Tならば 1 , その他(Cのみ)ならば0に符号化
df["fT"] = (df.f == "T").astype(int)

y = df.y
X = df[["x", "fT"]]
# 定数列(1)を作成
X = sm.add_constant(X)

# モデル生成と学習
model = sm.GLM(y, X, family=sm.families.Poisson())
result = model.fit()

# 結果出力
print(result.summary())

# 以下出力結果
                 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                      y   No. Observations:                  100
Model:                            GLM   Df Residuals:                       97
Model Family:                 Poisson   Df Model:                            2
Link Function:                    log   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -235.29
Date:                Sun, 21 Apr 2019   Deviance:                       84.808
Time:                        16:36:24   Pearson chi2:                     83.8
No. Iterations:                     4   Covariance Type:             nonrobust
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          1.2631      0.370      3.417      0.001       0.539       1.988
x              0.0801      0.037      2.162      0.031       0.007       0.153
fT            -0.0320      0.074     -0.430      0.667      -0.178       0.114
==============================================================================

本のP58に載っているのと全く同じ係数を得ることができました。

ついでに実際のデータとこのモデルから予測される種子数をプロットして可視化してみます。


import matplotlib.pyplot as plt
import numpy as np

# プロット用のデータ作成
xx = np.linspace(7, 13, 101)
yy0 = np.exp(result.params["const"] + result.params["x"]*xx)
yy1 = np.exp(result.params["const"] + result.params["x"]*xx + result.params["fT"])

# 可視化
fig = plt.figure(figsize=(8, 7))
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel("個体のサイズ")
ax.set_ylabel("種子の個数")
ax.scatter(X[X.fT == 0]["x"], y[X.fT == 0].values, marker="x", label="肥料無し")
ax.scatter(X[X.fT == 1]["x"], y[X.fT == 1].values, alpha=0.7, label="肥料有り")
ax.plot(xx, yy0, linestyle="--",  label="肥料無し")
ax.plot(xx, yy1, label="肥料有り")
plt.legend()
plt.show()

出力結果がこちら。

pandasで散布図行列を書く

何かしらのデータについて調べる時、各データごとの関係をまとめて可視化するために
散布図行列(matrix of scatter plots)を描きたくなることがあります。

matplotlibでゴリゴリ実装するのは手間なので、seabornが良く使われていると思うのですが、
実はpandasでも作成できるのでその紹介です。
(seabornを使えばいいじゃん、というのはごもっともなのですが、
なぜか僕の自宅環境ではwarningが発生し、その原因調査中なのであまり使いたくないのです。)

ドキュメントはこちら。
pandas.plotting.scatter_matrix

早速ですが、irisでやってみましょう。


import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris

iris = load_iris()
df = pd.DataFrame(
        iris.data,
        columns=iris.feature_names,
    )
pd.plotting.scatter_matrix(df, figsize=(15, 15), c=iris.target)
plt.show()

そして表示される結果がこちらです。

オプションで、 引数cにクラス情報を渡して色もつけました。
このコーディング量でかける図としては十分なクオリティではないでしょうか。

機械学習のプロダクト適用におけるリアルを語る会に参加しました

先日、ABEJAさんのオフィスで開催された、機械学習のプロダクト適用におけるリアルを語る会に参加しました。

AI/MLシステム開発の難しさ

発表者はディー・エヌ・エーの鈴木翔太さん。

企画からリリースまでの各段階において大変だったことを話されていました。
自分も経験があることが多くて全面的に納得です。

企画段階では、チームで何を作るのかのゴールの認識合わせが重要。
精度何%くらいを目指すのかと言った話もしておく必要があります。
また、全てを機械学習でやろうとしないことも大切ですね。

組織については、ML/AI開発チームにもソフトウェアエンジニアを入れて連携しているという話以前に、
MLチーム自体を作れてることを羨ましいと思いました。
(自分の職場には機械学習専門チームはなく、アナリストタスクと掛け持ちなので。)
懇親会で伺った話によるとやはりkaggler枠などの影響もあり、リファラルで採用が進んでいるそうです。

データ収集やアノテーション、データ管理が大変なのはどこでも語られている通り。
インフラをコスパ良く作るのも大変です。

また、誰が本番用コードを書くのかの話題も取り上げられました。
データサイエンティストが自分で書くのか、エンジニアに依頼するのか。
どちらも一長一短あります。
データサイエンティストのコードが汚いのは世界共通なんだろうか。

品質担保の話も取り上げられていたのは珍しいと感じました。
ただ、自分としてもここは悩んだ覚えがあります。
推論を含むシステムがエラーなく動くことをテストするのはいいとして、
予測である以上は結果が外れのこともあるわけで、
ここのデータに対して予測を外すのが正解みたいなテストを書くのも変で、
何がベストなのかよくわからなくなります。
そして最後はリリース後の評価やコストの話。

Deep Learningでリアルタイムに マーケット予測をしてみた

発表者はAlpaca Japanの梅澤慶介さん。

自分自身が個人トレーダーでもあるので密かに楽しみにしていたマーケット予測の事例。
ただ、やはり苦労されていたのはチームで開発するにあたって発生するあれこれのようでした。
メンバーの開発環境を再現できないとか、渡されたnotebookが動かないとか。
うちもここのメンバーのローカル環境揃ってないので良くわかります。
(自分一人に限ってもローカルPCで動いてたコードが開発用に建ててもらったインスタンスで動かないことがある。)
それをコンテナなどの技術で解決していったそうです。
Notebook禁止は衝撃。

機械学習のプラットフォームを自前で構築したり、
コードを共有してライブラリ化したりと、チームでの開発が順調に進んでいるようです。
ModelPackageみたいなの、僕も作りたいですね。

懇親会でも機械学習以外の部分含めて色々と質問ができ、非常に有意義でした。

智を集約しツラみを乗り越えたリピート推定の開発現場

発表者はABEJAの中川裕太さん。

題材はいくつかのミートアップですでに伺っているABEJA Insight for Retailのリピーター推定の話でした。
ただ、これまでと少し違った視点でのつらみの話があり興味深かったです。

「無限のビジネス可能性の前に立ち尽くす」ってのは新しい視点ですね。
たしかに、これができればいろんなことができそうっていうワクワクはありますが、
なんでもできるはなんにもできないと同じ、っていうのはあると思います。
(自分もかなり身に覚えがある)
顧客が片付けたいジョブに集中することが大事です。

研究と開発の時間軸が違って大変という話は、以前も伺った通りで
自分も直面するところなので良くわかります。
また、作ってる側から新しい手法が発見されるというも辛いところ。

コストと精度のトレードオフもありますね。
ここは、小売業向けのビジネスとしてのシビアな側面もありそうです。
(金融や医療、公的機関や大企業などお金があるところ向けでも無縁ではないと思いますが。)

今回はそれぞれの解決策もあり参考になりました。

マイクロサービスのDAGにして共通部分を使い回すなど、
実際、やってみるとかなり大変なのですがよく実現されたと思います。
スケーリングなどによるコスト最適化も、
アイデアは単純ですがおそらく実現は結構難しかったのではないかと。
精度の最適化の方にまだまだ課題も残っているようなので、今後に期待です。

まとめ

つらみの共有的な話を聞くことが増えていますが、
他社の人たちもそれぞれ苦戦しながら課題に立ち向かわれているのを知ると励みになりますね。
(普段は下手すると自分だけが効率悪く苦労しているような気がしてくるので。)

若干人数少なめな回だったのもあって、
懇親会でも発表者の方含めてゆっくり話ができて非常に有意義な時間を過ごせました。

機械学習の活用における課題は、
データサイエンティストだけで解決できるものだけではないと思うので、
この辺の知見はきちんと職場で共有して他チームの協力を得て取り組んでいきたいです。

Mac book に R環境構築

普段のデータ分析に使う環境はSQL, python, Tableau, Excelなどで完結させていたのですが、
業務上の都合でRも使うことになったので改めて学習を始めました。

Rを使うこと自他は初めてではなくpython使い始める前に何度か勉強したのですが、
どうにも好きになれなかったんで最近遠ざけていたのですがいよいよそうは言ってられなくなりましたね。

とりあえず、私物のMacにはインストールしていなかったので環境構築から始めます。
まずはRのインストールから。
配布サイトからダウンロードする方法もありますが、
homebrewを使う方法を採用しました。
いくつかのページを見ると、何かリポジトリを tapしないといけないように書いてあるところもあるのですが、
何もせずに brew install コマンドだけで入りました。


~$ brew install r
==> Installing dependencies for r: mpfr, gcc, libpng, openblas, pcre and readline
#(略) 依存パッケージとr本体がインストールされる

# 実行と終了
~$ r

R version 3.5.3 (2019-03-11) -- "Great Truth"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin18.2.0 (64-bit)

R は、自由なソフトウェアであり、「完全に無保証」です。
一定の条件に従えば、自由にこれを再配布することができます。
配布条件の詳細に関しては、'license()' あるいは 'licence()' と入力してください。

R は多くの貢献者による共同プロジェクトです。
詳しくは 'contributors()' と入力してください。
また、R や R のパッケージを出版物で引用する際の形式については
'citation()' と入力してください。

'demo()' と入力すればデモをみることができます。
'help()' とすればオンラインヘルプが出ます。
'help.start()' で HTML ブラウザによるヘルプがみられます。
'q()' と入力すれば R を終了します。

> q()

続いて、開発環境であるRStudioもインストールします。
これも配布元からダウロードする手もあるのですが、
homebrewで入れました。

これはアプリケーションフォルダに入れてGUIで使いたいので、brew cask を使います。


~$ brew cask install rstudio
==> Tapping homebrew/cask
Cloning into '/usr/local/Homebrew/Library/Taps/homebrew/homebrew-cask'...
remote: Enumerating objects: 4110, done.
remote: Counting objects: 100% (4110/4110), done.
remote: Compressing objects: 100% (4103/4103), done.
remote: Total 4110 (delta 24), reused 406 (delta 5), pack-reused 0
Receiving objects: 100% (4110/4110), 1.32 MiB | 187.00 KiB/s, done.
Resolving deltas: 100% (24/24), done.
Tapped 1 command and 4002 casks (4,116 files, 4.2MB).
==> Caveats
rstudio depends on R. The R Project provides official binaries:

  brew cask install r

Alternatively, the Homebrew-compiled version of R omits the GUI app:

  brew install r

==> Satisfying dependencies
==> Downloading https://download1.rstudio.org/desktop/macos/RStudio-1.2.1335.dmg
######################################################################## 100.0%
==> Verifying SHA-256 checksum for Cask 'rstudio'.
==> Installing Cask rstudio
==> Moving App 'RStudio.app' to '/Applications/RStudio.app'.
🍺  rstudio was successfully installed!

これで完了です。
念の為ですが、homebrew何も問題が起きてないことも見ておきました。


~$ brew doctor
Your system is ready to brew.

和分過程とARIMA過程

前回に引き続き今回も言葉の定義です。出典はいつもの沖本本から。

和分過程
d-1階差分をとった系列は非定常過程であるが、d階差分をとった系列が定常過程に従う過程は、
d次和分過程、もしくはI(d)過程と呼ばれます。
また、I(0)過程は定常過程で定義されます。

ARIMA過程
d階差分をとった系列が定常かつ反転可能なARMA(p,q)過程に従う過程は、
次数(p,d,q)の自己回帰和分移動平均過程、もしくはARIMA(p,d,q)過程と呼ばれます。

和分過程やARIMA過程の和分次数dは、AR特性方程式におけるz=1という解の個数に等しいことが知られているそうです。

前回の記事で定義した単位根過程は I(1)過程になり、
単位根過程の差分系列が定常で反転可能なARMA(p,q)過程の時、
それはARIMA(p,1,q)過程になります。

単位根過程の例として代表的なものにランダムウォークがあります。
言葉だけなら市場系のデータについて学んでいるとしょっちゅう出てきますね。

過程$y_t$が、
$$
y_t = \delta + y_{t-1} + \varepsilon_t , \ \ \ \varepsilon_t\sim iid(0, \sigma^2)
$$
と表現される時、$y_t$はランダムウォーク(random walk)と呼ばれます。
ただし、 $y_0=0$とします。また、 $\delta$はドリフト率と呼ばれるそうです。

定義式をよく見ると、AR(1)過程の形をしており、
AR特性方程式の解は$1$なので、これが定常でないこともわかります。
また、$\Delta y_t = y_t – y_{t-1} = \delta + \varepsilon_t$
であり、これは定常過程なので、ランダムウォークが単位根過程であることがわかります。

差分系列と単位根過程

久々に沖本先生の経済・ファイナンスデータの計量時系列分析に戻って、
時系列分析の用語を二つ紹介します。

まずは差分系列(5ページ)。
原系列$y_t$に対して、1点離れたデータとの差を取ることで生成できる系列、$y_t-y_{t-1}$を、
差分系列、または階差系列と呼び、$\Delta y_t$と書きます。

続いて単位根過程(105ページ)。
原系列$y_t$が非定常過程であり、差分系列$\Delta y_t=y_t-y_{t-1}$が定常過程である時、
過程は単位根過程である(unit root process)と言われます。

名前の由来は、単位根過程を誤差項が定常過程であるAR過程を用いて表現した時に、
AR特性方程式が$z=1$という解を1つ持つことらかきているそうです。

単位根過程の定義においては、差分系列が定常過程であることしか定義されておらず、
これがAR過程やARMA過程であることなどは要求されていないので注意です。

(お恥ずかしい話ですが、自分はデータサイエンティストに転職するずいぶん前に時系列分析を雑に勉強して、
差分とったらAR過程になるのが単位根過程だと誤って理解していたことがあります。)

このブログの通常の流れでは、ここから単位根過程のサンプルを一個構築して
プロットした図を出すのですが、実はすでに登場しているのでそちらにリンクしておきます。
参考:自己回帰過程のサンプル
この記事で構築したAR(3)過程の3個目の例が単位根過程です。

statsmodelsで重回帰分析

前回に引き続きpythonで重回帰分析を行い、回帰係数と切片を求める方法の紹介です。
今回はstatsmodelsを使います。

ドキュメントはこちら。
statsmodels.regression.linear_model.OLS

データは同じものを使い、結果が一致することを確認したいので
保存してたものを読み込みます。


import numpy as np
import statsmodels.api as sm

# データの読み込み
npzfile = np.load("sample_data.npz")
X = npzfile["X"]
y = npzfile["y"]

つぎに、回帰分析に入りますが、statsmodelsで回帰分析する場合には、一点注意があります。
それは、上で読み込んだ Xとyをそのままモデルに食わせると、切片(定数項)を0として結果を返してくることです。
それを回避するために、 X に対して、全ての値が1の列を追加する必要があります。
それは、次の関数を使って実行しています。
statsmodels.tools.tools.add_constant


# 配列Xに1だけからなる列を追加する。
X0 = sm.add_constant(X)
# 回帰分析実行
st_model = sm.OLS(y, X0)
results = st_model.fit()

# 切片と回帰係数
print(results.params)
# [ 4.94346432  2.00044153 -2.99255801  3.98231315  0.01708309]

無事にscikit-learnで行った時と全く同じ結果を得ることができましたね。
これだけだと、add_constantをやらなくてよかった分scikit-learnの圧勝に見えますが、
statsmodels を使うと一つメリットがあります。
それは、各係数が0でないことの検定ができることです。
resultsの sumary()って関数を実行してみましょう。


print(results.summary())

# 以下出力です。
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.999
Model:                            OLS   Adj. R-squared:                  0.999
Method:                 Least Squares   F-statistic:                 2.136e+04
Date:                Mon, 15 Apr 2019   Prob (F-statistic):          2.24e-139
Time:                        00:59:22   Log-Likelihood:                -143.94
No. Observations:                 100   AIC:                             297.9
Df Residuals:                      95   BIC:                             310.9
Df Model:                           4                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          4.9435      0.108     45.726      0.000       4.729       5.158
x1             2.0004      0.020    101.534      0.000       1.961       2.040
x2            -2.9926      0.018   -166.158      0.000      -3.028      -2.957
x3             3.9823      0.017    231.575      0.000       3.948       4.016
x4             0.0171      0.019      0.914      0.363      -0.020       0.054
==============================================================================
Omnibus:                        2.170   Durbin-Watson:                   1.994
Prob(Omnibus):                  0.338   Jarque-Bera (JB):                1.862
Skew:                          -0.208   Prob(JB):                        0.394
Kurtosis:                       2.477   Cond. No.                         7.12
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

出力の真ん中にある通り、定数項および回帰係数のt値と、p値が出力されています。
そして、x4 は P値が 0.05より大きいので、x4の係数が0であるという帰無仮説を棄却できていないことがわかります。
もともと、x4の回帰係数は0としてサンプルデータを作っていたので、これは良い結果です。
そのほか、AICなどの値も出してくれています。

scikit-learnでこれらの情報を収集するのは手間なので、
目的によってはstatsmodelsが便利な場面もあります。

scikit-learnで重回帰分析

今回と次回でpythonで重回帰分析を実行する方法を二つ紹介します。
今回はscikit-learnのLinearRegressionを使う方法です。

ドキュメントはこちら。
sklearn.linear_model.LinearRegression

最初に検証用のダミーデータを作ります。
$x_{i,j}$を -10 ~ 10の一様分布からサンプリングし、次の式で$y_i$を作ります。
$x_{i,3}$の係数が0になっていることから分かる通り、$x_{i,3}$は$y_i$には何の関係もない値です。
また、ノイズとして正規分布に従う乱数加えておきます。
$$
y_i = 5 + 2x_{i,0} -3x_{i,1} + 4x_{i,2} + 0x_{i,3} + \varepsilon_i,\\
\varepsilon_i \sim N(0,1) , \ \ \ i = 0,1,\cdots, 99
$$

サンプルデータを作って保存するコードがこちら。


import numpy as np
X = np.random.uniform(-10, 10, size=(100, 4))
y = 5 + X@[2, -3, 4, 0] + np.random.normal(0, 1, 100)
np.savez("sample_data.npz", X=X, y=y)

早速、回帰分析して回帰係数と定数項ついでに決定係数を求めてみましょう。
(回帰分析の目的が予測モデルを作ることであれば、データを訓練用と評価用に分けるべきなのですが、
今回は回帰係数を求める方法の紹介が主目的なので分けていません。)


import numpy as np
from sklearn.linear_model import LinearRegression

# データの読み込み
npzfile = np.load("sample_data.npz")
X = npzfile["X"]
y = npzfile["y"]

# モデルのインスタンス生成
model = LinearRegression()
#学習
model.fit(X, y)
# LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

# 回帰係数
print(model.coef_)
# [ 2.00044153 -2.99255801  3.98231315  0.01708309]

# 切片
print(model.intercept_)
# 4.943464324307898

# 決定係数 R^2
model.score(X, y)
# 0.9988893054003364

各コードの後ろにつけているコメントが僕の環境で実行した時の結果です。
回帰係数も切片もそれぞれほぼ正解に近い値が算出されていますね。

注:サンプルデータを乱数で生成しているので、データ生成からやり直せば結果は変わります。

scikit-learnを使うと、非常に手軽に回帰分析ができることがわかりました。
次回はstatsmodelsで同じことをやってみます。

numpyの配列をファイルに保存する

日常の実際で必要になったことはないのですが、
numpyに配列(ndarray)をファイルに保存する機能があるのでその紹介です。

(実用上、ファイルに保存したくなるようなデータは pandasのデータフレームに変換したあと、
to_csvやto_excelを使ってcsvかエクセルにしてしまうことが多いです。)

ドキュメントはこの辺り
numpy.save
numpy.load
numpy.savez

簡単に言ってしまえば、
numpy.save(ファイル名, 保存したい配列) で保存して、
numpy.load(ファイル名) で読み込める。
numpy.savez(ファイル名, (名前=)保存したい配列, (名前=)保存したい配列, …) で、
1ファイルに複数配列保存できる、と言った使い方になります。(名前=)は省略可能。

実際にやってみましょう。
まず検証用にダミーデータを作ります。


import numpy as np

# 保存のテスト用データを作成
data0 = np.arange(12).reshape(3, 4)
data1 = data0 ** 2
print(data0)
'''
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
'''
print(data1)
'''
[[  0   1   4   9]
 [ 16  25  36  49]
 [ 64  81 100 121]]
'''

まずは、save関数で配列を一つ保存し、それを別の変数に読み込んでみましょう。
特に難しいことはなく非常に直感的な使い方で上手く動きます。


# .npyフォーマットで保存
np.save("data0.npy", data0)

# 読み込み
data0_load = np.load("data0.npy")
print(data0_load)
'''
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
'''

続いて、savez関数で複数配列を1ファイルに保存し、復元してみます。
こちらは、loadした時に返されるオブジェクトが保存した配列そのものではないので注意です。
まずは名前をつけずに保存するパターン。


# savez を使うと複数の配列を.npzフォーマットの1ファイルにまとめて保存できる
np.savez("data0_1.npz", ary0, ary1)

# .npzファイルを読み込むと、配列ではなく、NpzFile オブジェクトが返される。
npzfile = np.load("data0_1.npz")
# .filesプロパティを見ると、中に二つの配列を持っていることがわかる
print(npzfile.files)
# ['arr_0', 'arr_1']

# それぞれの配列の取り出し
print(npzfile["arr_0"])
'''
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
'''
print(npzfile["arr_1"])
'''
[[  0   1   4   9]
 [ 16  25  36  49]
 [ 64  81 100 121]]
'''

また、名前をつけて保存すると次のようになります。


# 名前(例ではxとy)をつけて保存することも可能
np.savez("named_data0_1.npz", x=ary0, y=ary1)
# .npzファイルを読み込むと、配列ではなく、NpzFile オブジェクトが返される。
npzfile2 = np.load("named_data0_1.npz")
# .filesプロパティを見ると、中に二つの配列を持っていることがわかる
print(npzfile2.files)
# ['x', 'y']

# 保存時の名前で取り出すことが可能
print(npzfile2["x"])
'''
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
'''
print(npzfile2["y"])
'''
[[  0   1   4   9]
 [ 16  25  36  49]
 [ 64  81 100 121]]
'''

勝手にarr_0とか命名されるよりも、自分で名前をつけておいた方が混乱がなさそうですね。