2020-10-08

APIキー、OAuthクライアントID、サービスアカウントキーの違い:Google APIs

GoogleのサービスをAPI経由で利用する際、大きく3つの認証情報が登場します。「APIキー」、「OAuth2.0クライアントID」、「サービスアカウントキー」です。本投稿ではこれらの用途をざっくりと整理し、取得方法、Pythonでの使い方の違いを説明します。

f:id:messefor:20201006223438p:plain

まとめ
基本事項のおさらい
1. APIキー
- 取得方法
- Pythonでの利用例
2.OAuth2.0クライントID
- 取得方法
  - 1.同意画面の設定
  - 2.OAuthクライアントIDの作成
- Pythonでの利用例
3. サービスアカウントキー
- 取得方法
- Pythonでの利用例

まとめ

Google APIsを使うには、下表の3ついずれかの認証情報が必要です
APIでどういうデータに、どういう形でアクセスするかに依って使うべき認証情報が異なります
いずれの認証情報もGoogle Cloud Consoleの APIとサービス > 認証情報から作成することができます
- Googleアカウントが無い方は、まずGoogleアカウントを作成し、Google Cloud Consoleでプロジェクト作成が必要です

#	必要となるもの	アクセス可能なデータ	誰としてアクセスするか	具体例
1	APIキー	一般公開データ	匿名ユーザ	YouTubeにある動画をアプリケーション経由で検索する
2	OAuth2.0クライアトID	一般公開データ/ユーザーデータ	ユーザアカウント	あるユーザ（エンドユーザ）の代わりにユーザのGoogleドライブにアプリケーションを経由でデータを保存
3	サービスアカウントキー	一般公開データ/ユーザーデータ	サービスアカウント	共同作業メンバのGoogleカレンダー情報にアプリケーションを経由してアクセスする

基本事項のおさらい

認証情報について理解する上で知っておいたほうが良い基本事項についてまとめます。必要ない方は飛ばしてください。

APIの性質

API（Application Programming Interface）はサービスを外部のプログラムやアプリケーションから利用できるようにした窓口です。例えば、Google MapはGoogleのサービスの一つで、普通私たちユーザはブラウザやGoogle Map Appを自分で操作してGoogle Mapを利用します。一方で自身が開発するアプリの中でGoogle Map機能の一部を利用したい場合は、アプリの中のプログラムを介してGoogle Map機能を操作する必要があります。ここで利用されるのがAPIというプログラムから呼び出すインターフェースです。

一般にAPIはプログラムやアプリケーションから利用されます。

認証情報の役割

Google は許可されたアプリケーションからしかAPI利用をさせてくれません。そのため我々API利用者はGoogleアカウントを作成し、Google Cloud Console内でAPIを利用するアプリケーションの登録を行う必要があります。

つまりどの開発者が登録した、どのアプリケーションがAPIを叩いているのかをGoolge様にお知らせするのがAPIキーやOAuth2.0クライアトIDなど認証情報の役割の一つです。

またAPI利用は使用量の無料枠上限を超えると課金が発生します。API利用がどのアプリでどの程度発生しているかをGoogle様にお知らせするのも認証情報の役割です。課金管理はプロジェクト単位で行われるため認証情報はGoogleアカウント内のプロジェクトと紐付ける必要があります。これが認証情報を作成する際にプロジェクトが必要な理由です。

3つの認証情報は想定される利用シーンが異なる

APIキー、OAuth2.0クライアントID、サービスアカウントキーの3種類のいずれかがあればAPIを使うことが可能です。ただし三者は用途が異なり、用途にあった認証情報を選択するのが自然です。たとえば一般公開データにしかアクセスしないのに、OAuthクライアントIDを使う必要はありません。用途に応じた必要最低限の認証情報を使います。

以上を踏まえて３つの認証情報をそれぞれ見ていきましょう。

1. APIキー

まずはAPIキーについてです。Google Cloud ドキュメント：API キーの使用には次のように説明されています。

API キーは、プリンシパルなしでアプリケーションを識別する暗号化された単純な文字列です。一般公開データに匿名でアクセスする場合に便利で、割り当てや課金のために API リクエストをプロジェクトに関連付けるために使用されます。

YouTubeでの検索機能のように一般公開データにのみアクセスできれば十分な場合、利用者にユーザアカウントでログインしてもらうような必要はありません。このようなケースではAPIキーを使ってAPIを利用します。APIキーは暗号化された単純な文字列で、作成も取り扱いも比較的簡単です。

取得方法

事前にGoogle Cloud アカウント作成とプロジェクト作成が必要です。まだな方は以下のサイトなどを参考にして、アカウントとプロジェクト作成を行ってください。

Google アカウントを作成する（Gmail アドレスも同時に作成） | MAGELLAN BLOCKS

プロジェクト作成が完了したら次の手順でAPIキーを作成します。

1. プロジェクトを選択

Google Cloud Consoleにログインして、APIを利用するプロジェクトを選択します。

f:id:messefor:20201005221836p:plain:w250

APIとサービス > 認証情報

f:id:messefor:20201005221902p:plain:w250

f:id:messefor:20201005221928p:plain:w250

2. APIキーの作成

「認証情報の作成」をクリック

f:id:messefor:20201005221959p:plain:w250

「APIキー」を選択

f:id:messefor:20201005222022p:plain:w250

APIキーが表示されたボックスの右側にあるコピーマークをクリックし、閉じる

f:id:messefor:20201005222040p:plain:w250

3. キーで利用するAPIを制限する（オプショナル）

作成直後の状態ではすべてのAPIの利用が可能となっているので、利用したいAPIのみに制限します

作成したAPIキーを選択

f:id:messefor:20201005222059p:plain:w250

APIの制限 > キーを制限

f:id:messefor:20201005222115p:plain:w250

Pythonでの利用例

例として、PythonからAPIキーを使ってYouTubeのsearchAPIを使うサンプルコードを掲載します。　下記コードのYOUTUBE_API_SERVICE_NAMEやYOUTUBE_API_VERSIONを変えるといろいろなAPIを使えるようになります。

from googleapiclient.discovery import build

# 利用するAPIサービス
YOUTUBE_API_SERVICE_NAME = 'youtube'
YOUTUBE_API_VERSION = 'v3'

# APIキーを指定
YOUTUBE_API_KEY = 'Your-API-KEY-here'

# API のビルドと初期化
youtube = build(YOUTUBE_API_SERVICE_NAME, YOUTUBE_API_VERSION,
            developerKey=YOUTUBE_API_KEY)

# YouTubeでキーワード検索（視聴数の多い順）
q = 'Overjoyed'
response = youtube.search().list(part='snippet',
                                q=q,
                                order='viewCount',
                                type='video',).execute()

response

>>>  # 出力

{'etag': 'QvSya4g44W8aUkApKibFr-o5JqA',
 'items': [{'etag': 'Jw9SPInmYcoJg21LlKvbdDhk1IU',
            'id': {'kind': 'youtube#video', 'videoId': '_a1LogyX9Uw'},
            'kind': 'youtube#searchResult',
            'snippet': {'channelId': 'UCOo4Oc-PYSrzEeappVyoM-Q',
                        'channelTitle': 'beechoppers',
                        'description': 'please enjoy the music! music and '
                                       'lyrics by stevie wonder clip design by '
                                       'R. Norton.',
                        'liveBroadcastContent': 'none',
...

2.OAuth2.0クライントID

アプリ経由で利用者のGoogleドライブにファイルを保存をする場合など、アプリケーションがエンドユーザーに代わって Google Cloud APIs にアクセスします。このケースで利用するのがOAuth2.0クライアントIDです。OAuth2.0プロトコルを利用することで、アプリケーションが特定ユーザとしてサービスを利用します。その際に許可する権限の範囲（スコープ）も制限します。

f:id:messefor:20201005230723p:plain

OAuth2.0クライントIDは、OAuth2.0クライアントクレデンシャル（資格情報）というJSONファイルとしてダウンロードし利用することが多いと思います。OAuth2.0クライアントクレデンシャルには、OAuth2.0認証を利用するアプリの情報（クライアントID、クライアントシークレット）が含まれます。クレデンシャルはGoogle様にアプリ開発者がOAuth認証を行う資格があることを示します。

取得方法

OAuth認証は作成する同意画面の設定やアプリの種類など状況によって適切に設定する必要がありますが、本投稿ではテストとして最短で利用する手順を書いてみます。

1.同意画面の設定

OAuth認証では、認証時ユーザにログインをしてもらったり、許可してもらう権限の範囲を示す必要があるため、認証画面に表示する情報は重要です。以下はconnpassにTwitterアカウントでログインする際のOAuth認証画面です。画面にはconnpassのロゴなど表示され、どのアプリがアカウント情報を利用しようとしているかひと目で分かります。

f:id:messefor:20201005230918p:plain:w400

このユーザに同意を求める部分の設定を「同意画面の設定」から行っていきます。今回はテストなのでなるべく手軽に行うことを目指します。

APIキー同様、まずGoogle Cloud Consoleにログインして、APIを利用するプロジェクトを選択します。

f:id:messefor:20201005224733p:plain:w250

APIとサービス > OAuth同意画面を選択し、適当なアプリ名とメールアドレスを入力します

f:id:messefor:20201005225025p:plain:w250

UserTypeを聞かれるので、今回は「外部」を選びます。「内部」は使ったことありませんが、G Suiteの組織内で利用するときに使うようです。

f:id:messefor:20201005225113p:plain:w250

省略できる箇所は省略し、試すのに必要最低限な入力をして「保存して次へ」で進めます。入力が必要なのはメールアドレスのみです。実際にアプリを開発する際はスコープなど慎重に設定しましょう。

f:id:messefor:20201005225301p:plain:w250

f:id:messefor:20201005225358p:plain:w250

f:id:messefor:20201005225404p:plain:w250

f:id:messefor:20201005225409p:plain:w250

2.OAuthクライアントIDの作成

「＋認証情報を作成」から「OAuthクライアントID」を選択します f:id:messefor:20201005225647p:plain:w250

アプリケーションの種類はここでは、デスクトップを選択しています f:id:messefor:20201005230153p:plain:w250

これでOAuthクライアントクレデンシャルが作成されました f:id:messefor:20201005225857p:plain:w250

「OK」を押して、作成されたIDの右にある↓矢印からJSONをダウンロードします。

f:id:messefor:20201005230357p:plain:w250

これでOAuthクライアントIDを含む、OAuthクライアントクレデンシャルが発行されました。

Pythonでの利用例

例として、PythonからOAuthクレデンシャルを使ってYouTubeのplaylistAPIを使うサンプルコードを掲載します。以下はOAuth認証用の関数を定義をしています。この関数を使ってOAuth認証をし、その後ユーザデータにアクセスします。今回はYouTubeに自分用の再生リストを作成したいと思います。実行のコードを下部にあります。

'''OAuth認証用の関数'''

import os
import pickle

import numpy as np
import pandas as pd

from googleapiclient.discovery import build
from google_auth_oauthlib.flow import InstalledAppFlow
from google.auth.transport.requests import Request


def get_credentials(client_secret_file, scopes,
                    token_storage_pkl='token.pickle'):
    '''google_auth_oauthlibを利用してOAuth2認証

        下記URLのコードをほぼそのまま利用。Apache 2.0
        https://developers.google.com/drive/api/v3/quickstart/python#step_1_turn_on_the_api_name
    '''
    creds = None
    # token.pickleファイルにユーザのアクセス情報とトークンが保存される
    # ファイルは初回の認証フローで自動的に作成される
    if os.path.exists(token_storage_pkl):
        with open(token_storage_pkl, 'rb') as token:
            creds = pickle.load(token)

    # 有効なクレデンシャルがなければ、ユーザーにログインしてもらう
    if not creds or not creds.valid:
        if creds and creds.expired and creds.refresh_token:
            creds.refresh(Request())
        else:
            flow = InstalledAppFlow.from_client_secrets_file(
                client_secret_file, scopes=scopes)
            creds = flow.run_local_server(port=0)

        # クレデンシャルを保存（次回以降の認証のため）
        with open(token_storage_pkl, 'wb') as token:
            pickle.dump(creds, token)

    return creds

以下のコードのCLIENT_SECRET_FILEに先程ダウンロードしたJSONを指定してください。これを実行するとまずOAuth認証のためにブラウザが立ち上がり利用者にGoogleアカウントでのログインを求めます。認証が終わると「再生リストテスト」という名で空の再生リストが作成されます。

'''OAuth認証と再生リストの作成'''
# 利用するAPIサービス
YOUTUBE_API_SERVICE_NAME = 'youtube'
YOUTUBE_API_VERSION = 'v3'

# OAuthのスコープとクレデンシャルファイル
YOUTUBE_READ_WRITE_SCOPE = 'https://www.googleapis.com/auth/youtube'
CLIENT_SECRET_FILE = 'client_secret.json'

# OAuth認証：クレデンシャルを作成
creds = get_credentials(
                    client_secret_file=CLIENT_SECRET_FILE,
                    scopes=YOUTUBE_READ_WRITE_SCOPE
                    )

# API のビルドと初期化
youtube_auth = build(YOUTUBE_API_SERVICE_NAME, YOUTUBE_API_VERSION,
                    credentials=creds)

# Add Playlist
# This code creates a new, private playlist in the authorized user's channel.
playlists_insert_response = youtube_auth.playlists().insert(
  part="snippet,status",
  body=dict(
    snippet=dict(
      title='再生リストテスト',
      description='APIで作成したプレイリスト'
    ),
    status=dict(
      privacyStatus="private"
    )
  )
).execute()

3. サービスアカウントキー

サービスアカウントはアカウントの一種です。アカウントなのでメールアドレスが割り当てられますが、アプリケーションが利用するアカウントなのでユーザアカウントとは異なり、ログインなどはできません。利用者がこのアプリケーション用のアカウントとデータを共有することで、プログラムからユーザデータへのアクセスが可能になります。APIはこのサービスアカウントの代わってサービスの操作を行います。

用途としてはたとえば、自分のGoogleカレンダーをサービスアカウントと共有し、アプリケーションからAPI経由で閲覧や編集するなどがあります。

　サービスアカウントキーでの認証は、一般にサービスアカウントクレデンシャルというJSONファイルを使って行われます。サービスアカウントクレデンシャルにはアカウントのメールアドレスや、鍵情報、クライントIDなど認証に必要な情報が含まれています。これまでの認証と違ってアカウント認証もこの資格情報で行います。

取得方法

1. プロジェクトを選択

Google Cloud Consoleにログインして、APIを利用するプロジェクトを選択します。

f:id:messefor:20201005221836p:plain:w250

APIとサービス > 認証情報

f:id:messefor:20201005221902p:plain:w250

f:id:messefor:20201005221928p:plain:w250

2. サービスアカウントキーの作成

「+認証情報を作成」から「サービスアカウント」を選択 f:id:messefor:20201007213935p:plain:w250

サービスアカウント名を入力

f:id:messefor:20201007213940p:plain:w250

続行を選んで、完了します。

f:id:messefor:20201007213944p:plain:w250 f:id:messefor:20201007213950p:plain:w250

3. 鍵ペアの作成

作成したサービスアカウントを選んで f:id:messefor:20201007215021p:plain:w250

鍵を追加で、「新しい鍵を追加」してJSONを選びます。ついでに、上部のメールアドレスもコピーしておきましょう。

f:id:messefor:20201007220531p:plain:w250

f:id:messefor:20201007214957p:plain:w250

保存した.jsonファイルを認証に利用します。

Pythonでの利用例

例として、Pythonからサービスアカウントキーを用いてGoogle スプレッドシートのセルの値を取得してみます。まずClould Consoleの各プロジェクトでGoogle Sheet APIを有効にしておきましょう。

APIとサービス > ライブラリ

f:id:messefor:20201007215809p:plain:w250

検索窓でsheetを入力して、Google Sheets APIを選ぶ

f:id:messefor:20201007215843p:plain:w250

Google Sheets APIを有効にします

f:id:messefor:20201007215444p:plain:w250

次の準備として、Google スプレッドシートで新規にスプレッドシートを作成し、作成したシートをサービスアカウントと共有（編集権限を与える）しておく必要があります。

何か入力して、右上の「share」を押下

f:id:messefor:20201007220147p:plain:w250

さきほどメモったサービスアカウントのメールアドレスを貼り付けて「OK」を押して、シート共有します。

f:id:messefor:20201007220754p:plain:w250

さらにシートのURLから以下のsheet_idをメモっておきましょう https://docs.google.com/spreadsheets/d/<sheet_id>

実行するコードこちらです。これでシートの値が取得できます。

from googleapiclient.discovery import build
from google.oauth2 import service_account

CREDENTIAL_FILE = 'credential_json_path'
credentials =\
  service_account.Credentials.from_service_account_file(CREDENTIAL_FILE)

SAMPLE_SPREADSHEET_ID = '1CL0pg_LDpASaHQqagDxBDKonmD64NRS6CQnpAbVDjhs'
SAMPLE_RANGE_NAME = 'Sheet1!A1:B4'

# Sheets APIビルド
service = build('sheets', 'v4', credentials=credentials)

# セルの値を取得
sheet = service.spreadsheets()
result = sheet.values().get(spreadsheetId=SAMPLE_SPREADSHEET_ID,
                            range=SAMPLE_RANGE_NAME).execute()
values = result.get('values')


>>>  # 出力

[['1', 'あ'], ['2', 'い'], ['3', 'う']]

APIで権限まわりは混乱しがちですね。ちょっと長くなりましたが整理してみました。日々精進

2020-10-04

YouTube API: 再生時間が条件に合う動画を選んで再生リストを作成（2/2）

Python YouTube Data API YouTube B-life

前回はB-lifeのYouTubeチャンネルを題材に、チャンネル中にある動画情報一覧および再生時間を取得し、再生時間が10分以上15分以下の条件に収まる動画のみに絞り込みを行いました。本投稿では、これら動画を使って再生リストを作成します。

前回掲載した全ステップの中のでいうと本投稿はステップ3を説明します。

~~特定のチャンネル（B-life）の中にある動画ID一覧を取得する~~ （前回済）
~~動画の再生時間を取得し、再生時間で絞り込む~~ （前回済）
自分のYouTubeアカウント上に該当する動画の再生リストを作成する
1. 準備
2. OAuth認証とAPIのビルド
3. 再生リストの作成と動画の追加

messefor.hatenablog.com

1.準備

APIの認証情報について

前回searchAPIを使う際はAPIキーさえあれば良かったですが、API経由でYouTube再生リストを作成するにはOAuth 2.0 クライアントクレデンシャルが必要です。これは再生リストの作成・編集の際、ユーザデータへのアクセスが必要なためです。ユーザデータへのアクセスはOAuth2.0認証で行われますが、OAuth認証のためにOAuthクライアントクレデンシャルが必要になります。OAuth2.0認証はユーザーデータへ安全にアクセスできるよう、許可する権限範囲など細かに指定できる仕組みを提供します。

今回具体的に必要なものは、クライアントクレデンシャル情報が書き込まれた.jsonファイルです。Google API Consoleで作成してダウンロードしておきます。

APIキーやOAuthクライアントキーについての基本的な内容は以下を参考にしてください。 messefor.hatenablog.com

認証用ライブラリのインストール

認証部分はGoogleのドキュメント：quickstart.pyにあるコードの丸パクリですが、このサンプルコードがgoogle_auth_oauthlibを使っているのでインストールしておきます。

pip google-auth-oauthlib

2. OAuth認証とAPIのビルド

OAuth認証の処理部分は以下のget_credentials関数で行います。引数token_storage_pklで指定された場所に有効なトークンがなければ、ユーザにログイン（OAuth認証処理）をしてもらい、データへのアクセスを許可してもらいます。ユーザの許可により有効なトークンが得られれば、次回以降のアクセスのためにそれをtoken_storage_pklにpickleとして保存します。

def get_credentials(client_secret_file, scopes,
                    token_storage_pkl='token.pickle'):
    '''google_auth_oauthlibを利用してOAuth2認証

        下記URLのコードをほぼそのまま利用。Apache 2.0
        https://developers.google.com/drive/api/v3/quickstart/python#step_1_turn_on_the_api_name
    '''
    creds = None
    # token.pickleファイルにユーザのアクセス情報とトークンが保存される
    # ファイルは初回の認証フローで自動的に作成される
    if os.path.exists(token_storage_pkl):
        with open(token_storage_pkl, 'rb') as token:
            creds = pickle.load(token)
            
    # 有効なクレデンシャルがなければ、ユーザーにログインしてもらう
    if not creds or not creds.valid:
        if creds and creds.expired and creds.refresh_token:
            creds.refresh(Request())
        else:
            flow = InstalledAppFlow.from_client_secrets_file(
                client_secret_file, scopes=scopes)
            creds = flow.run_local_server(port=0)
            
        # クレデンシャルを保存（次回以降の認証のため）
        with open(token_storage_pkl, 'wb') as token:
            pickle.dump(creds, token)

    return creds

実行部分のコードは以下です。YOUTUBE_READ_WRITE_SCOPEがOAuth認証で使われるスコープになります。今回はYouTubeアカウントの管理権限（全読み書き）を範囲にしていますが、OAuth 2.0 Scopes for Google APIsを見て、適宜必要な最小限のスコープを絞ってください。

'''OAuth認証とAPIのビルド実行'''

# 利用するAPIサービス
YOUTUBE_API_SERVICE_NAME = 'youtube'
YOUTUBE_API_VERSION = 'v3'

# OAuthのスコープとクレデンシャルファイル
YOUTUBE_READ_WRITE_SCOPE = 'https://www.googleapis.com/auth/youtube'
CLIENT_SECRET_FILE = 'config/client_secret.json'

# OAuth認証：クレデンシャルを作成
creds = get_credentials(
                    client_secret_file=CLIENT_SECRET_FILE,
                    scopes=YOUTUBE_READ_WRITE_SCOPE
                    )

# API のビルドと初期化
youtube_auth = build(YOUTUBE_API_SERVICE_NAME, YOUTUBE_API_VERSION,
                    credentials=creds)


>>> 
# 出力されるメッセージ
Please visit this URL to authorize this application: https://accounts.google.com/o/oauth2/auth?response_type=code&client_id=hoge.apps.googleusercontent.com&redirect_uri=http%3A%2F%2Flocalhost%3A54887%2F&scope=https%3A%2F%2Fwww.googleapis.com%2Fauth%2Fyoutube&state=hoge&access_type=offline

このコードを実行すると、token_storage_pklで指定された場所に有効なトークンがない場合、ブラウザが立ち上がって以下のようなGoolge ユーザログイン画面が表示されます。（ブラウザが立ち上がらない場合は、出力されるメッセージに含まれるURLをクリックして、OAuth認証を行ってください）

f:id:messefor:20201004124432p:plain:w300

さらに警告画面が表示されますが、実験なので無視して詳細 >（安全でないページ）に移動と遷移すると

f:id:messefor:20201004124527p:plain:w300

YouTubeアカウントの管理権限を許可するかどうか聞かれますので、許可を選択します。

f:id:messefor:20201004124551p:plain:w300

再度確認されるので、許可を押します。

f:id:messefor:20201004124657p:plain:w300

これでユーザ側のOAuth認証許可の処理が完了になります。API経由で再生リストを作成するための権限が与えられました。次は実際に再生リストを作成して、動画を追加します。

3. 再生リストの作成と動画の追加

新規再生リストの作成

まず新規再生リストを追加します。タイトルと説明、再生リストを公開にするか限定するかを指定しています。

'''新規再生リストの作成'''

title = 'B-life Test 10分〜15分'
description = '再生時間が10分〜15分の動画の再生リスト'
privacy_status = 'public'  # 'private'

# 新規再生リストを追加
# https://developers.google.com/youtube/v3/docs/playlists/insert
playlists_insert_response = youtube_auth.playlists().insert(
  part="snippet, status",
  body=dict(
    snippet=dict(
      title=title,
      description=description
    ),
    status=dict(
      privacyStatus=privacy_status
    )
  )
).execute()

これを実行して、YouTubeサイトに行って自分のライブラリを確認するとB-life Test 10分〜15分という再生リストが作成されています。

f:id:messefor:20201004124807p:plain:w300

まだ動画は一本も含まれていないので、クリックしても「この再生リストには動画がありません」と表示されます。

f:id:messefor:20201004124741p:plain:w300

再生リストへの動画の追加

追加する動画は、前回は作成した動画IDリストとします。df_video_playlistデータフレームのidカラムに動画IDが格納されています。また、追加先となる再生リストのIDは新規再生リストの作成時のレスポンスとして取得できていますので、それを利用します。

'''再生リストへの動画の追加 '''

# 動画ID
videoids = df_video_playlist['id'].values

# example videodis
# 前回の続きとして実行しない場合は以下のコメントを外す
# videoids = ['nZSe3ZZUSJw', 'EYqPH91q5nE', 'Zw_osQctQKs', 'kBLCmz8pODg',
#             'IBt5l9Um_rY', 'WWxKW8ncPe0', 'CmWJJWKLeKg', 'e9HR7-3I3MA',]

# プレイリストに動画を追加
# https://stackoverflow.com/questions/20650415/insert-video-into-a-playlist-with-youtube-api-v3/22190766

playlistid = playlists_insert_response['id']  # 作成した再生リストのIDを取得

# 動画IDをループ
for videoid in videoids:
  
    resourceid = dict(kind='youtube#video',
                      videoId=videoid)

    response = youtube_auth.playlistItems().insert(
                part='snippet',
                  body=dict(
                    snippet=dict(
                      playlistId=playlistid,
                      resourceId=resourceid
                      )
                   )
            ).execute()

先程の再生リストを見ると10分〜15分の動画が追加されているのが分かります。Mariko先生がバッチリ映えております。コード全体はGithubにアップしています。

f:id:messefor:20201004124852p:plain:w300

YouTube Data APIで再生リストを作成し動画を追加しました。動画コンテンツの整理は時間だけでなく、自然言語処理を使ったカテゴリ分類や音声や動画を用いたラベル付けなど、色々拡張ができそうです。日々精進

2020-10-04

YouTube API: 再生時間が条件に合う動画を選んで再生リストを作成（1/2）

Python B-life Google API YouTube YouTube Data API

YouTubeで再生時間の条件で動画の再生リストを作りたくなることありませんか？今回YouTube Data APIを使って再生時間が10分以上、15分以下の動画のみを含む再生リスト作成するPythonコードを書いたので共有します。長くなりそうなので、2回に分けて投稿します。

messefor.hatenablog.com

背景：YouTubeのフィルタが微妙

運動不足解消のため、毎朝妻と一緒にB-lifeというYouTubeのヨガレッスンやっています。B-lifeのチャンネルには質の高いレッスン動画が多くアップロードされていて嬉しいのですが、毎朝やっているとその日やる動画を選ぶのも結構迷います。このとき「今日は時間がないので5分のレッスンが良い」「今日は15分」というような状況が頻繁に発生するので、再生時間の条件で再生リストが作れると便利です。

ただ残念ながらYouTubeの検索フィルタでは細かい再生時間の絞り込みはできなさそうです。

f:id:messefor:20201004001021p:plain:w300

YouTube Data APIを使ったらできるだろうということでやってみました。

実現までのステップ

YouTube Data APIを使った今回の実現ステップをまとめると次の3つです。

特定のチャンネル（B-life）の中にある動画ID一覧を取得する
動画の再生時間を取得し、再生時間の条件に合う動画を選ぶ
自分のYouTubeアカウント上に該当する動画の再生リストを作成する

1と2に関しては、認証も必要なくAPIキーさえ用意すれば気軽に行えます。3にはOAuth認証が必要になります。本投稿では、第1回として上記1と2の手順を一つひとつ見ていきたいと思います。3は次の投稿で説明します。

APIキーやOAuthクライアントキーについては以下を参考にしてください。

messefor.hatenablog.com

0. 準備

YouTubeチャンネルのIDを確認

APIで動画やチャンネルを指定するためには、それらのIDが必要です。準備として対象とするYouTubeチャンネルのIDを確認します。チャンネルIDはURLから直接分かります。YouTubeチャンネルのURLはhttps://www.youtube.com/channel/<チャンネルID>となっているので、ブラウザでアクセスして確認します。

ちなみにB-lifeのチャンネルIDはUCd0pUnH7i5CM-Y8xRe7cZVgでした。

B-lifeのYouTubeチャンネルのホーム画面

赤枠がチャンネルIDです。 f:id:messefor:20201004001106p:plain:w400

なおB-lifeは初心者向けの動画も多く、Mariko先生のレクチャが大変分かりやすいのでおすすめです

ライブラリをインストール

pipでPythonのGoogle API Client、をインストールしておきましょう。あとpandasも使いますのでなければインストールしましょう。

pip install --upgrade google-api-python-client
pip install pandas

1. チャンネルにある動画ID一覧を取得

さっそくYouTubeチャンネルの中にある動画ID一覧を取得していきます。

APIのビルドと初期化

まずAPIの初期化です。 Google APIを使うためには、APIキーが必要です。Google Cloud Consoleから作成して、下記のYOUTUBE_API_KEY = 'your-api-key-here'を作成したAPIキーに書き換えてください

from googleapiclient.discovery import build

# 利用するAPIサービス
YOUTUBE_API_SERVICE_NAME = 'youtube'
YOUTUBE_API_VERSION = 'v3'

# APIキー
YOUTUBE_API_KEY = 'your-api-key-here'

# API のビルドと初期化
youtube = build(YOUTUBE_API_SERVICE_NAME, YOUTUBE_API_VERSION,
                developerKey=YOUTUBE_API_KEY)

動画情報一覧を取得

次にsearch()APIを使って動画ID一覧を取得するのですが、それを関数get_video_list_in_channel()で定義しています。 get_video_list_in_channel()はチャンネルIDを入力とし、動画情報を出力します。仕様としては、指定したチャンネルIDの中で公開日時が新しいものから50件ずつ最大2回リクエストされます。なので最大50 * m_re_cntの動画情報が取得されます。

Google APIでは1日あたりのAPI使用量（クォータ）が限られているので、沢山リクエストすると一発で上限を超えてしまいます。今回実験ということで引数max_req_cnt=2と小さめに指定しています。必要に応じて調整してください。

from datetime import datetime, timedelta
import pandas as pd

        
def get_video_list_in_channel(youtube, channel_id, max_req_cnt=2):
    '''特定のチャンネルの動画情報を取得し、必要な動画情報を返す

        公開時刻が新しい順に50ずつリクエスト
        デフォルトでは最大2リクエストで終了
    '''

    n_requested = 50

    earliest_publishedtime =\
        datetime.now(tz=timezone.utc).strftime('%Y-%m-%dT%H:%M:%SZ')

    req_cnt = 0
    result = []
    while True:
        response = youtube.search().list(part='snippet',
                                        channelId=channel_id,
                                        order='date',
                                        type='video',
                                        publishedBefore=earliest_publishedtime,
                                        maxResults=n_requested).execute()
        req_cnt += 1
        video_info = fetch_video_info(response)
        result.append(video_info)

        # 取得動画の最も遅い公開日時の1秒前以前を次の動画一覧の取得条件とする
        last_publishedtime = video_info['publishTime'].min()
        last_publishedtime_next =\
            datetime.strptime(last_publishedtime, '%Y-%m-%dT%H:%M:%SZ') - timedelta(seconds=1)

        earliest_publishedtime =\
            last_publishedtime_next.strftime('%Y-%m-%dT%H:%M:%SZ')

        if req_cnt > max_req_cnt:
            # リクエスト回数がmax_req_cntを超えたらループを抜ける
            print('Result count exceeded max count {}.'.format(max_req_cnt))
            break

        if len(response['items']) < n_requested:
            # リクエストした動画数より少ない数が返った場合はループを抜ける
            print('Number of results are less than {}.'.format(n_requested))
            break

    if len(result) > 1:
        df_video_list = pd.concat(result, axis=0).reset_index(drop=True)
    else:
        df_video_list = result[0]

    return df_video_list


def fetch_video_info(response, as_df=True):
    '''APIのレスポンスから必要な動画情報を抜き出す'''
    info_list = []
    for item in response['items']:
        info = {}
        info['title'] = item['snippet']['title']
        info['kind'] = item['id']['kind']
        info['videoId'] = item['id']['videoId']
        info['description'] = item['snippet']['description']
        info['publishTime'] = item['snippet']['publishTime']
        info['channelTitle'] = item['snippet']['channelTitle']
        info['thumbnails_url'] = item['snippet']['thumbnails']['default']['url']
        info_list.append(info)
    if as_df:
        return pd.DataFrame(info_list)
    else:
        info_list

実行部分のコードは以下です。df_video_listに取得した情報が格納されています。df_video_listのvideoId列が動画IDになります。動画ID一覧が取得できました。

# 動画一覧を取得したいチャンネルID
channel_id = 'UCd0pUnH7i5CM-Y8xRe7cZVg'

# 動画一覧を取得
df_video_list = get_video_list_in_channel(youtube, channel_id)

# 出力
df_video_list.head()

f:id:messefor:20201004001151p:plain

2. 動画の再生時間を取得し、フィルタをかける

上で取得した動画情報一覧に動画の再生時間は含まれていませんので、別途videoIdをキーにしてAPI経由で取得する必要があります。再生時間を取得した後はpandasで条件に合ったものを絞り込むだけです。

再生時間を取得

再生時間を取得するにはvideoAPIに動画IDを渡して詳細情報を取得します。ここでは、下表のような機能のラッパー関数をいくつか定義しています。

関数	機能概要
get_contents_detail(), get_contents_detail_core()	動画の詳細情報取得
get_duration(), get_basicinfo(),	レスポンスjsonのパース
pt2sec()	再生時間のフォーマット変換

1回のリクエストで50動画IDまでしか問い合わせできないみたいなので、get_contents_detail()では動画ID一覧を50件ずつに分割してから処理しています。

import re
import numpy as np

def get_contents_detail_core(youtube, videoids):
    '''動画の詳細情報を取得'''
    part = ['snippet', 'contentDetails']
    response = youtube.videos().list(part=part, id=videoids).execute()
    results = []
    for item in response['items']:
        info = get_basicinfo(item)
        info['duration'] = get_duration(item)
        results.append(info)
    return pd.DataFrame(results)


def get_contents_detail(youtube, videoids):
    '''必要に応じて50件ずつにIDを分割し、詳細情報を取得'''
    n_req_pre_once = 50

    # IDの数が多い場合は50件ずつ動画IDのリストを作成
    if len(videoids) > n_req_pre_once:
        videoids_list = np.array_split(videoids, len(videoids) // n_req_pre_once + 1)
    else:
        videoids_list = [videoids,]

    # 50件ずつ動画IDのリストを渡し、動画の詳細情報を取得
    details_list = []
    for vids in videoids_list:
        df_video_details_part =\
            get_contents_detail_core(youtube, vids.tolist())
        details_list.append(df_video_details_part)

    df_video_details =\
        pd.concat(details_list, axis=0).reset_index(drop=True)
    return df_video_details


def get_duration(item):
    '''動画時間を抜き出す（ISO表記を秒に変換）'''
    content_details = item['contentDetails']
    pt_time = content_details['duration']
    return pt2sec(pt_time)


def get_basicinfo(item):
    '''動画の基本情報の抜き出し'''
    basicinfo = dict(id=item['id'])
    # snippets
    keys = ('title', 'description', 'channelTitle')
    snippets = {k: item['snippet'][k] for k in keys}
    basicinfo.update(snippets)
    return basicinfo


def pt2sec(pt_time):
    '''ISO表記の動画時間を秒に変換 '''
    pttn_time = re.compile(r'PT(\d+H)?(\d+M)?(\d+S)?')
    keys = ['hours', 'minutes', 'seconds']
    m = pttn_time.search(pt_time)
    if m:
        kwargs = {k: 0 if v is None else int(v[:-1])
                    for k, v in zip(keys, m.groups())}
        return timedelta(**kwargs).total_seconds()
    else:
        msg = '{} is not valid ISO time format.'.format(pt_time)
        raise ValueError(msg)

再生時間取得の実行部分です。df_video_detailsが結果で、durationカラムに再生時間が入っています。再生時間の単位は秒です。

# 動画ID
videoids = df_video_list['videoId'].values

# 動画の詳細情報を取得
df_video_details = get_contents_detail(youtube, videoids)

df_video_details.head()

f:id:messefor:20201004001214p:plain

再生時間でフィルタをかける

ここではpandas.between()を使って10分以上15分以下の動画に絞り込みこみました。

# 再生時間が10分以上、15分以下の動画に絞り込む
lower_duration = 10 * 60  # 10分以上
upper_duration = 15 * 60  # 15分以下
is_matched = df_video_details['duration'].between(lower_duration, upper_duration)

df_video_playlist = df_video_details.loc[is_matched, :]

# 出力
df_video_playlist.head()

これで10分以上、15分以下の動画IDのリストができました。めでたし、めでした。次回はこれら動画を含む再生リストをAPI経由で生成したいと思います。コードはGithubに挙げています。

久々にGoogle APIを触りました。サービスも年々増えているみたいで、やりさえすれば色々できそうです。やりさえすれば。。。日々精進

2020-09-12

ファイルをバイト単位で比較：cmpコマンド

macOS Bash シェルスクリプト

2つのzipファイルを比較したい場合など、バイナリレベルでファイルを比較したいことがあります。こういう場合はシェルのcmpコマンドで比較すると良いです。

FILE1とFILE2を比較する場合の構文は以下になります。

cmp FILE1 FILE2

完全に一致していれば何も出力されず、不一致の場合は

file1.tar.gz file2.tar.gz differ: char 5, line 1

のように比較した各ファイルパスとどの部分が異なるのかを標準出力に出力してくれます。（ここではfile1.tar.gzとfile2.tar.gzを比較しています）

また差異の有無によって終了ステータスコードが異なり、差があった場合はexit status = 1、差がなかった場合は exit status = 0となります。差の有無によって分岐処理などしたい場合はexit statusで判断すれば良いかなと思います。

例：2つのディレクトリにある全ファイルをバイト単位で比較する

使用例として、2つのディレクトリにある同名のファイル同士を比較し、差があるときだけメッセージを出すbash スクリプトを書いてみます。差の有無の判断はexit statusで行っています。

#! /bin/bash
# 2つのディレクトリにある同じ名前のファイル同士を比較し、差があるときだけメッセージを出す

dir1=src1/  # １つ目のディレクトリ
dir2=src2/  # ２つ目のディレクトリ

fpaths=${dir1}*  # 第１ディレクトリの全ファイルパスを取得

# ファイルパスをループ
i=0
for f1path in $fpaths; do
  fnm=`basename $f1path`
  # echo $fnm
  f2path=${dir2}${fnm}
  cmt=`cmp $f1path $f2path`  # ２つのファイルを比較
  if [ $? -eq 1 ]; then
    echo $cmt  # 差があった場合は、メッセージを表示
    i=$(( $i + 1 ))  # 差があったファイル数をカウントしておく
  fi
done

# 差が１つもなかった場合はその旨を表示
if [ $i -eq 0 ]; then
  echo 'All file pairs are the same.'
fi

検証環境：

macOS Catalina

GNU bash, version 3.2.57(1)-release (x86_64-apple-darwin19) Copyright (C) 2007 Free Software Foundation, Inc.

コードはこちらです。

shellとかRとか、少し触らないとすぐに忘れてしまいます。栓が抜けた湯船のように、蛇口から入ってくる水より多くの水が抜けていきます。そのうちスッカラカンになりそうな。。。日々精進

2020-09-07

Softbank光が遅いのでWiFiルータを買ったら爆速になった話

光回線通信機器 SoftBank光 WiFiルータ

Softbank光に加入して数ヶ月が経ちました。レンタル品（光BBユニットEWMTA2.3）のWiFi/ルータ機能を使っていたのですが、どうも通信が安定せずにストレスフルな状態でした。そこでWiFiルータを新しく調達したところ大幅に快適になりました。私の環境でどのように改善したかまとめます。 f:id:messefor:20200907232634p:plain

今回分かったこと（結論）
改善前の状況
- 参考：なぜ光BBユニットを使っているか
- 改善前の速度
改善策と手順
- 購入したルータ
どの程度改善したか
まとめ

今回分かったこと（結論）

環境に依ると思いますが、光BBユニット使っている人は市販のWiFiルータ買うと通信ストレスが大幅に改善される気がします。

光BBユニットのWiFiルータは電波の飛びが弱く不安定
- 特にルータから遠い部屋では通信が遅い。5GHz帯が弱い印象
- 速度計測だけでは測れない不安定さがある
市販のWiFiルータを調達することで、大幅に通信速度が改善し、快適になる
バッファロー AirStation WSR-1800AX4/NWHは比較的安価で個人的にはおすすめでした

改善前の状況

5月にSoftbank光を導入し「これで爆速になる」とワクワクしていたのだが、実際に使い始めるとどうも様子がおかしい。

我が家ではルータ（AP：アクセスポイント）はリビングに置いていて、私の作業部屋からAPまで扉を2つ挟んで10mくらいの距離がある。作業部屋でリモート会議をしていると、かなり頻繁に音声が途切れる。平均すると20秒に1-2秒くらい聞こえないんじゃないかというひどい状況。新型コロナの影響で在宅勤務が増えている中、仕事にも支障が出てきてたので対策をすることにした。

参考：なぜ光BBユニットを使っているか

携帯がY!mobileなのでSoftbank光 + おうち割光セット が適用されると携帯の通信費も安くなる（端末を2台使っているので毎月1000円以上割引）。このおうち割の適用条件に「光BBユニットレンタル」が含まれているため、レンタルして使っていた。

光BBユニットレンタル

「Wi-Fiマルチパック」または「Wi-Fi地デジパック」

「ホワイト光電話」、「BBフォン」または「光電話（N）およびBBフォン」

つまり光BBユニットをレンタルすることでおうち割（1000円割引）を手に入れているのだ。特に光BBユニットを使いたいわけではない。そんなもんだろう。おうち割についてはこのサイトが分かりやすい

改善前の速度

作業部屋での速度計測サイトで10回ずつ計測すると、帯域・上り下りの平均はだいたい30Mbpsくらい。ただ5GHz帯が上り下りともに弱いのがグラフでも見て取れる。

f:id:messefor:20200907232802p:plain

※エラーバーは±1標準誤差

さらに計測すると、光BBユニットWiFiルータが置いてあるリビングでは140Mbps、有線では295Mbpsと十分な速度がでていることも分かった。さすがに有線だとリモート会議をしていても全く途切れない。音質も良いとのこと。以上からAPから離れた場所でのWiFiの速度低下が問題だろうと考え、そこを改善しようと考えた。

改善策と手順

WiFiの改善アプローチだが、以下2つを思いた。

WiFiルータ（アクセスポイント）を購入する
WiFi中継機を購入する

結構迷ったが以下の記事に背中を押され、まずはWiFiルータを購入することにした。設定手順もこの記事に全て書いてある。ありがたく参考にした。

ソフトバンク光に無線LANルーターを追加したら、爆速安定して人生が多少変わった

購入したルータ

バッファロー AirStation WSR-1800AX4-WHを購入。

www.buffalo.jp

値段は8,000円くらいまでで検討した。他の選定条件は以下

項目	条件	理由
価格	8000円くらいまで	家庭の事情
色	ホワイト	家庭の事情
発売時期	比較的新しい機種	電化製品は新しい方が良いと思っているから
機能	ビームフォーミング	改善前の電波が弱く改善したいから
	WiFi 6対応	今後使うかもしれないから
	IPv6	今後使うかもしれないから
	中継機能	今後使うかもしれないから

ルータはほとんどが黒のものが多く、色と値段だけでかなり絞り込まれる。また中継機能がついているルータならば、もし今回のルータ交換で速度が改善しなくても作業部屋に中継機として設置できると考えた。

どの程度改善したか

5GHz帯がめちゃくちゃ改善した。速度で見ても3倍程度出ている。2.4Hz帯は数値ではさほど変わりがないが、2.4GHzも使っていて安定感が大幅にアップした印象を受ける。どちらも通信がほぼ途切れなくなった。リモート会議でも音声が途切れることはほぼなくなり、ストレスから開放された。Webページ開くときのレスポンスも結構早くなった印象を受ける。

結局我が家の環境では、5GHz帯が遠くまで飛ぶことが重要だったと考えている。これがルータを買うことで改善されたとみた。

f:id:messefor:20200907233052p:plain

5GHz帯、有線の速度比較@作業部屋

	アクセスポイント	5GHz帯平均速度	リモート会議や動画通話
WiFi（改善前）	光BBユニットEWMTA2.3	27.2Mbps	頻繁に途切れる
WiFi（改善後）	WSR-1800AX4/NWH	84.3Mbps	問題ない（たまに途切れることがある）
有線	-	295.5Mbps	途切れない。音が良い（空気感が伝わる）

表を見ても分かるとおり、なんだかんだで有線は圧倒的に通信速度が早い。有線だと空気感が伝わるほど高音質とのことで、会話相手にも評判だ。なお私はMacBookProを使ってリモート会議など参加しているのですが、マイクの性能はそこそこ良いみたい。通信が良いと、マイクやスピーカも買おうかなという気持ちになってくる。。。

まとめ

光BBユニットのWiFiルータは電波の飛びが弱く不安定
- 特にルータから遠い部屋では通信が遅い。5GHz帯が弱い印象
- また速度計測だけでは測れない不安定さがある
Softbank光では市販のWiFiルータを調達することで、大幅に速度が改善し快適になる
- 5GHzが遠くても比較的強い
BUFFALO AirStation WSR-1800AX4/NWHは安価で個人的にはおすすめでした
けっきょく有線最強

在宅勤務では通信環境が精神的な負担をかなり左右します。改善してから随分ストレスフリーになりました。同じような症状の方は市販ルータの購入を検討してみてもいいと思います。

技術以外での日常の奮闘もどこかに出力しておきたいなと思って書きました。日常分からない事だらけですね。日々精進

2020-09-05

Stanでガンマ回帰（動かす編）

stan R データ分析ベイズ推定

回帰問題で目的変数が正の連続値をとる場合、ガンマ回帰は選択肢の一つになります。個人的に適用シーンは多いと思うのですが比較的情報が少ない気がします。本投稿ではトイデータとStanでサクッとモデルを動かしてみたいと思います。なおガンマ回帰はGLM（一般化線形モデル）の枠組みで推定できるためRのglm()でも簡単にfitできます。シンプルにモデリングしたいだけであれば、あえてStanでやる意味はないかもしれません

設定
トイデータ生成
Stanファイル
MCMC実行
パラメタ推定結果
推定したパラメタからのサンプリング
まとめ
参考文献

設定

今回の設定として目的変数 $\displaystyle{y_i}$ の期待値 $\displaystyle{\mu_i}$ が、説明変数 $\displaystyle{x_i}$ のべき乗の線形和 $\displaystyle{b_{0} + b_{1}x_{i}^ {a}}$ で表現できるとします。

$\displaystyle{ E[y_i] = \mu_{i} = b_{0} + b_{1}x_{i}^{a} }$

推定するパラメタは $\displaystyle{b_0, b_1, a}$ です。あえて $\displaystyle{a}$ を入れたのは、Stanなら冪乗の項の対数をとるなど工夫せずに直接推定できそうだと思ったからです。またGLMでいうリンク関数はidentity（恒等リンク関数）です。そのため $\displaystyle{\mu_i \leq 0}$ とならないように説明変数の範囲を注意しないといけません。

さて期待値 $\displaystyle{\mu_i}$ でガンマ分布に従う $\displaystyle{y_i}$ を考えたいのですが、ガンマ分布はshapeパラメタ: $\displaystyle{\alpha}$ とrateパラメタ: $\displaystyle{\beta}$ によって次のように決まります。（ $\displaystyle{y, \alpha, \beta}$ は正の実数）

$\displaystyle{ \operatorname{Gamma}(y|\alpha, \beta) = \frac{\beta^{\alpha}}{\Gamma(\alpha)}y^{\alpha-1}\operatorname{exp}^{-\beta y} }$

われわれに馴染み深い期待値 $\displaystyle{\mu}$ と分散 $\displaystyle{\phi}$ は、 $\displaystyle{\alpha, \beta}$ で次のように表現できます。正規分布と違って $\displaystyle{\mu}$ と $\displaystyle{\phi}$ は互いに強く依存しています。

$\displaystyle{ \begin{align} \mu &= \frac{\alpha}{\beta} \tag{1} \\ \phi &= \frac{\alpha}{\beta^2} \tag{2} \end{align} }$

ガンマ分布では $\displaystyle{\alpha}$ と $\displaystyle{\beta}$ を決める必要があるため、 $\displaystyle{(1)(2)}$ より $\displaystyle{\alpha}$ と $\displaystyle{ \beta}$ を $\displaystyle{\mu}$ と $\displaystyle{\phi}$ で表します。

$\displaystyle{ \begin{align} \alpha &= \frac{\mu^2}{\phi} \\ \beta &= \frac{\mu}{\phi} \\ \end{align} }$

今回の設定では $\displaystyle{\mu}$ を説明変数の関数で表します。また $\displaystyle{\phi}$ はGLMの枠組みでは一定と仮定するのが一般的なようです。今回も同様の仮定を置きます。以上が今回のデータ生成過程の設定になります。モデリングではこれらの関係性と観測データをStanに直接指示すれば推定できそうです。まずはトイデータを生成しましょう。

トイデータ生成

Rを使って学習データを生成します。 $\displaystyle{a = 0.3, b_0 = 1.0, b_1 = 1.5, \phi = 0.1}$ としています。

set.seed(1234)

N <- 150

# 推定したいパラメタ
a <- 0.3
b0 <- 1.0
b1 <- 1.5
phi <- 0.1

# 一様分布からxを生成
x <- runif(n=N, 0.1, 20)

# 期待値と線形予測子
mu <- b0 + b1 * x ^ a

# パラメタの変換
shape <- mu ^ 2 / phi
rate <- mu / phi

# ガンマ分布に従う乱数生成
y <- rgamma(n = length(x), shape = shape, rate = rate)

生成したデータをプロットします。 $\displaystyle{x}$ と $\displaystyle{y}$ の散布図です。

library(tidyverse)

# x vs. y の散布図プロット
g <-list(x = x, y = y, mu = mu) %>% data.frame() %>% ggplot()
g <- g + geom_point(aes(x=x, y=y))
g <- g + geom_line(aes(x=x, y=mu, color='mu'))
g <- g + ggtitle('Toy data $mu = 0.2 + 1.5 * x ^ 0.3'))
g

f:id:messefor:20200905230011p:plain

Stanファイル

先程のデータ生成過程をそのままStanファイルに記載するだけです。これをexample01.stanとして保存します。

data {
  int N;  // 観測数
  vector[N] x;  // 説明変数x
  vector[N] y; // 目的変数y
}

parameters {
    // 推定したいパラメタ
  real<lower=0> a;
  real<lower=0> b0;
  real<lower=0> b1;
  real<lower=0> phi;  // variance of gamma
}

transformed parameters {
  vector<lower=0>[N] alpha;  // ガンマ分布のshapeパラメタ
  vector<lower=0>[N] beta;  // ガンマ分布のratioパラメタ
  vector<lower=0>[N] mu;  // 期待値

  for (n in 1:N) {
    mu[n] = b0 + b1 * pow(x[n], a);
  }

  for (n in 1:N) {
    alpha[n] = pow(mu[n], 2.) / phi;
  }
  
  beta = mu / phi;
}

model {
  y ~ gamma(alpha, beta); // likelihood
}

generated quantities {
  vector<lower=0>[N] y_new;
  for (n in 1:N) {
    y_new[n] = gamma_rng(alpha[n], beta[n]);
  }
}

MCMC実行

rstanを使って事後分布からのランダムサンプリングを実行します。

library(rstan)

options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

data <- list(N = N, x = x, y = y)
fit <- stan(file = 'example01.stan', data = data, seed = 1234, iter = 2000)

divergent transitionsの警告が出ますが、警告の遷移数が少ないのでここでは無視します。

Warning message: “There were 2 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup” Warning message: “Examine the pairs() plot to diagnose sampling problems

警告の詳細は以下を参考にしてください。 messefor.hatenablog.com

収束も問題なさそうです。

f:id:messefor:20200905230148p:plain

パラメタ推定結果

どのパラメタも真の値 $\displaystyle{a = 0.3, b_0 = 1.0, b_1 = 1.5, \phi = 0.1}$ を分布の~~25-75%~~20-80%の間に含んでいるので、大きく外していないようです。

	mean	se_mean	sd	2.5%	25%	50%	75%	97.5%	n_eff	Rhat
a	0.2643425	0.0018063278	0.04469031	0.20230515	0.22890738	0.2563850	0.2925427	0.3675151	612.1158	1.0053146
b0	0.6953202	0.0173214550	0.43094516	0.02972518	0.33524315	0.6733556	1.0266480	1.5374601	618.9781	1.0035314
b1	1.8591395	0.0164183020	0.40759706	1.06712580	1.54369900	1.8705854	2.1966878	2.5121268	616.3191	1.0038491
phi	0.1051311	0.0003396671	0.01246250	0.08323707	0.09635969	0.1042309	0.1122532	0.1334601	1346.1805	0.9992741

f:id:messefor:20200905230237p:plain

推定したパラメタからのサンプリング

推定したパラメタからのサンプリングデータを使って、学習データと乖離していないか確認します。

# 推定パラメタからのサンプリングを抽出
is.pred <- str_detect(rownames(result.summay), 'y_new.')
data <- data.frame(result.summay[is.pred,])

colnames(data) <-
  c('mean', 'se_mean', 'sd', 'p2.5', 'p25', 'p50',
      'p75', 'p97.5', 'n_eff', 'Rhat')
data$x <- x
data$y <- y
data$mu <- mu

# 図示
g <- ggplot(data=data)
g <- g + geom_line(aes(x=x, y=mean, color='post_mean')) +
      geom_line(aes(x=x, y=mu, color='mu')) +
      geom_point(aes(x=x, y=y))

g <- g + geom_ribbon(aes(x=x, ymin=p25,ymax=p75), fill="blue", alpha=0.2) +
      geom_ribbon(aes(x=x, ymin=p2.5,ymax=p97.5), fill="blue", alpha=0.2)
g <- g + ggtitle('True and Predicted') + labs(y='y')
g

水色の線が事後分布の期待値です。また濃い青が50%予測区間、薄い青が95%予測区間です。大きな乖離はなさそうですし、推定できているとしましょう。

f:id:messefor:20200905230330p:plain

コード全体はここに置いてあります。Rのコード、Stanのコード

まとめ

Stanを使ってガンマ回帰を行いました。正規分布を仮定した線形回帰での推定は、目的変数の分布の対称性が仮定されていることや負の値が出てきてしまうため、今回のようなケースではガンマ回帰の方がより自然な選択肢になります。GLMでも推定は行えますが、不確実性やモデルを柔軟に構築したい場合はStanを使うのも手かもしれません。

統計モデリングでは確率分布の知識が要求されますね。日々精進

参考文献

2020-08-29

stan: Divergent transitions after warmupの警告が出る場合

R stan MCMC ベイズ推定

stanでサンプリングを行っていると次のような警告がでる場合があります。

1: There were 15 divergent transitions after warmup.
2: Examine the pairs() plot to diagnose sampling problems

おっかないので、ドキュメントを読んで整理してみました。部分的にしか理解してない状態ですので、間違いがありましたら教えてもらえると嬉しいです。

なんの警告なのか

MCMCでは推定するパラメタを表す擬似粒子を考えて、事後分布が作り出す（位置エネルギー）場を粒子がどういうふうに動くか、その軌跡をシミュレートします。この軌跡がある一定の基準から遠くに大きく外れて行くことをdivergence（発散？逸脱？）と呼んでいます。divergenceが大きくなると、このシミュレーションは信頼できないため、devegent trasitionsの警告がでます。

どいういうときに発生するか

シミュレーションは、1ステップごと逐次的に探索をしながら行われるのですが、事後分布がぐにゃぐにゃにカーブしているケースでは、ステップ幅が事後分布のカーブに対して十分に小さな値でないと、分布の形状をキャプチャできず、本来たどるべき軌跡から逸脱してしまいます。一方でステップ幅が小さすぎると、探索が少しずつしか進まず非効率になったり、途中で（Uターンする前に）停止してしまう可能性があります。最急降下法と同じ感じですね。場所によってカーブの大きさが異なるような事後分布では、ある場所では非効率な探索ステップになり、ある場所では逸脱が大きくなるようなステップになってしまいます。これが警告の主な原因とのこと。

どうすればいいか

シンプルな対処としてadapt_deltaパラメタを大きくする と書いてあります。adapt_deltaはデフォルトで0.8に設定されていますので、stan()のcontrol引数でlist(adapt_delta = 0.99)のように0.8より大きい値を指定します。この値を大きくするとステップ幅を強制的に小さくできます。

# `adapt_delta`を0.8より大きくして、ステップ幅を小さくする
fit <- stan(..., control = list(adapt_delta = 0.99))

結果的にサンプリングに時間がかかるようになる可能性がありますが、ウォームアップ後にdivergenceが発生して結果の信頼性が落ちるよりはマシだろうという考えです。

もう少し詳細な対処法

警告文に出ているようにpairs()を使って結果を図示して診断することで、もう少し詳しく診断できます。

# 結果のペアプロット
pairs(fit)

f:id:messefor:20200829231632p:plain

見にくいかもしれませんが、グラフ中の小さな赤い点がdevegent trasitionsですが、この赤い点が1. 対角の下半分（下三角）のグラフに出ているか、2.上半分（上三角）に出ているか、で性質や対処法が若干変わります。上の例では下三角にdevegent trasitionsが出ています。

1.上三角にでている場合は、単純にadapt_deltaを大きくすれば良いですが、2.下三角に出ている場合adapt_deltaを大きくするだけで改善するケースは限られています。devegent trasitionsの数が小さい場合は改善する可能性があるので、とりあえずにadapt_deltaを大きくしましょう。それでも改善しない場合は、事後分布の形がシンプルになるように、モデルの表現を変更する必要があります。特に階層モデルを使う場合にこのケースが発生しやすいとのこと。

対処法をまとめると、

`devegent trasitions`（赤い点）が出る位置	対処
1.上三角部分にあるグラフに出ている	`adapt_delta` を大きくする
2.下三角部分にあるグラフに出ている	- `adapt_delta` を大きくする(改善するケースは限定的) - 事後分布の形が単純になるようにモデルを変更する

事後分布の形がシンプルになるように、モデルの表現を変更する というのはかなり難しそうに思えますが、別機会に詳細をまとめてみたいと思います。

物理出身の方は、MCMCやHMCをすんなり理解できるみたいですが、物理分かってない私みたいな人間には目が回りそうです。日々精進