- 投稿日：2019-10-20T23:43:54+09:00
AWS Toolkit for Visual Studio Codeを使ってみた（C#編）
AWS Toolkit for Visual Studio Codeとは
AWS Toolkit for Visual Studio CodeはAWSがオープンソースで開発を行っているVisual Studio Code（以下VS Code）用のプラグインで、ソースコードはGitHubのaws/aws-toolkit-vscode で公開され開発が進められています。 ライセンスは Apache License 2.0で、2019/10/20 時点の最新バージョンは1.1.0です。
AWS Toolkit for Visual Studio Codeの便利な機能
AWS Toolkit for Visual Studio Codeは、VS Codeのプラグインとして下記の便利な機能を提供します。
- AWS Serverless Application Model (以下 SAM) を利用したサーバーレスアプリケーションの開発
- AWS Lambdaの開発言語として、Node.js、Python、.NET Coreに対応
- SAMテンプレートを含むAWS Lambdaファンクションプロジェクトの生成
- ユニットテスト用コードの提供
- VS Code内でローカル実行やデバッグ実行
- VS Code内からSAMテンプレートによるAWSへのデプロイ
利用するための準備
前提条件
AWS Toolkit for Visual Studio Codeをインストールするために必要な前提条件は以下の通りです。
- アマゾンウェブサービスのアカウント
- サポートされている以下のいずれかのオペレーティングシステム
- Windows
- Linux
- macOS
- 使用する言語のSDK
- Visual Studio Code バージョン 1.31.1 以降
また使用する言語に応じて下記のVS Codeプラグインをインストールします。
AWSサーバーレスアプリケーション開発を行うには以下もインストールします。
筆者は下記の環境で動作を確認しました。
- macOS Mojave バージョン 10.14.6
- Visual Studio Code バージョン 1.39.2
- AWS Toolkit for Visual Studio バージョン 1.1.0
- C# for Visual Studio Code (powered by OmniSharp) 1.21.5
- .NET Core 3.0.100
- aws-cli/1.16.156 Python/3.6.1 Darwin/18.7.0 botocore/1.12.146
- SAM CLI, version 0.22.0
- Docker version 19.03.2, build 6a30dfc
AWS Toolkit for Visual Studio Codeのインストール
インストールは、シンプルにプラグインをインストールするだけです。
AWSで利用するプロファイルの設定
AWS Toolkit for Visual Studio Codeをインストールしたら、AWSのリソースにアクセスするために認証情報の設定を行います。コマンドパレットから [AWS: Create Credentials Profile]を選択します。
初期プロファイル名、利用するIAM Userの認証情報（アクセスキーIDとシークレットアクセスキー）を入力します。すでに設定されている場合は、~/.aws/credentialsファイルと ~/.aws/configファイルが開きますので、内容を確認します。アクセスキーIDとシークレットアクセスキーの取得の方法は、AWSアクセスキーの取得を参照してください。
AWS Toolkit for Visual Studio Codeの利用
AWS Toolkit for Visual Studio Codeのインストールとプロファイルの設定が完了したら、サーバーレスアプリケーションを開始することができます。
AWSサーバーレスアプリケーションの作成
サーバーレスアプリケーションの開発を始めるために、コマンドパレットから[AWS: Create New SAM Application]を選択します。
作成するサーバーレスアプリケーションの言語を選択します。今回はC#編ということで、dotnetcore2.1を選択します。
ワークスペースのフォルダを選択します。
アプリケーション名を入力します。
サーバーレスアプリケーションプロジェクトが作成され、template.ymlが表示されます。Amazon API Gatewayから呼び出されるAWS LambdaをAWS上に生成するSAMテンプレートが作成されていることがわかります。
ローカル環境でデバッグモードの実行
サーバーレスアプリケーションプロジェクトには、AWS Lambdaファンクションのサンプルコードが予め作成されています。srcフォルダ内のFunction.csを開きます。
左側のアクティビティーバーからAWSアイコンをクリックするとCodeLensが表示されます。デバッグ実行で実行を中断したい行番号の左側をクリックしてブレークポイントを設定します。AWS LambdaファンクションのエントリポイントであるFunctionFunctionHandlerメソッドのCodeLensで[Debug Loccaly]をクリックしてデバッグ実行を開始します。
デバッグツールバーとデバッグコンソールが表示され、.NET Coreアプリケーションがデバッグモードで実行開始されます。
デバッグツールバーの続行（F5）をクリックするとブレイクポイントを設定した箇所で実行が停止することを確認できます。またローカル変数の値が確認できることもわかります。
API Gatewayイベントソースリクエストからの情報の受け取り
Amazon API GatewayをイベントソースとするAWS Lambdaファンクションでイベントソースからの値を受け取るためにAWS Lambda for .NET Coreでは、Amazon.Lambda.APIGatewayEvents.APIGatewayProxyRequestクラスを提供しています。このクラスはパブリックプロパティを集めただけのクラスですが、AWS API Gatewayから渡されるJSONのメンバーの値が自動的にこのクラスのインスタンスにセットされて渡されるようになっています。この動作をデバッグモードの実行でエミュレートするための機能が提供されています。
CodeLensの[Configure]をクリックすると.aws/templates.jsonが開きます。
このevnetプロパティに値をセットするとデバッグ実行にその値を受け取ることができます。サーバーレスアプリケーションプロジェクト内の[アプリケーション]/eents/event.jsonにAWS API Gatewayから渡されるjsonのサンプルがあるので、これをコピーして試してみます。
再度、CodeLensの[Debug Locally]をクリックしてデバッグモードで実行して見るとAPIGatewayProxyRequestクラスのインスタンスであるapigProxyEventに値がセットされて渡されていることが、変数ウィンドウの値で確認できます。
AWSへのデプロイ
AWS Toolkit for Visual Studio Codeが提供する機能を利用して、VS CodeからサーバーレスアプリケーションをAWSへデプロイすることができます。デプロイする前に、あらかじめAmazon S3にデプロイ用のバケットを作成しておきます。
VS Codeのコマンドパレットから[AWS: Deploy SAM Application]を選択します。
SAMテンプレートを選択します。
次にデプロイするリージョンを選択します。
あらかじめ作成しておいたAmazon S3のバケット名を入力します。バケットは1つ前で選択したリージョンに作成している必要があります。
デプロイスタックのための名前を入力します。
以上の手順でVS Codeからサーバーレスアプリケーションをデプロイすることができます。
最後に
AWS Toolkit for Visual Studio Codeを利用することで使い慣れたVS Codeを使用してAWS Lambdaの開発を効率よく行うことができます。今回はC#（.NET Core）で開発を行う例を紹介しましたが、Node.jsやPythonでも同様に利用することができます。この投稿が、VS Codeでサーバーレスアプリケーション開発したいと考えている開発者の方の参考になれば幸いです。
免責
本投稿は、個人の意見で、所属する企業や団体は関係ありません。
また掲載しているサンプルプログラム等の動作に関しても一切保証しておりません。
- 投稿日：2019-10-20T15:27:34+09:00
Unity玉転がしチュートリアル 3-1.収集するオブジェクトの作成
この記事の対象者
- Unity入門したい人
- 最初の一歩が踏み出せない人
OSとか環境とか
- Windows 10 Pro
- macOS Mojave
- Unity 2019.2.8f1
- Rider 2019.2.2
補足
- 公式動画にて利用しているのはMacなので、Windowsユーザーはある程度脳内変換して見る事
- 筆者はWindows、Macの両方の環境で確認。Ubuntuとかでは検証してない。
- 基本Unityは英語メニューで利用
- 間違いがあったらツッコミ大歓迎
公式
https://learn.unity.com/tutorial/collecting-scoring-and-building-the-game#5c7f8529edbc2a002053b788
オブジェクトを作成
プレイヤーが獲得するアイテムを作成
GameObject > 3D Object > Cubeで追加→座標リセットを行う追加すると既存のプレイヤーのオブジェクトが操作に邪魔なので、非アクティブにする
非表示にする方法はInspectorにて名前の左側にあるチェックボックのON/OFFで切り替わる
追加したCubeが埋まっているので、PositionのYを0.5上に移動してあげる
メラミンスポンジみたいな形なったら成功
オブジェクトの外観変化
アイテムっぽくないので見た目を変更する
追加したアイテムの
・ScaleをXYZ全て0.5に
・TransformのRotationのXYZを全て45に
にして、小さく地面から浮いた状態に変更するこれだけでもアイテム感は確かにあるが、アイテムが動いていれば更にアイテムとわからせることが出来るので動きを追加
Cubeを回転させるスクリプトを追加
Pick Upオブジェクトを選択してInspectorからAdd ComponentでNewScriptを選択してスクリプト追加
出来上がったスクリプトはこちらRotator.csusing System.Collections; using System.Collections.Generic; using UnityEngine; public class Rotator : MonoBehaviour { // Update is called once per frame void Update() { transform.Rotate(new Vector3(15,30,45) * Time.deltaTime); } }これでPlayを押すと、アイテムがしっかり回転していることがわかる
Pick Up オブジェクトのプレハブ化
プレハブ＝GameObjectの設計図
プレハブを更新すると全てのオブジェクトが更新されるプレハブを格納するフォルダを作成する
フォルダ名は「Prefabs」とする
HierarchyのPick UpからProjectのPrefabsの中にドラッグするゲームオブジェクトの整理
整理用のオブジェクト（Pick Ups）を作って、その中を整理します。
Pick UPオブジェクトをPick Upsの子に設定する
このままPick Upオブジェクトを移動すると、斜め（今の角度準拠）に動いてしまう
行いたいことは地面にオブジェクトを動かす事なので、エディタのモードを変更するLocal→Globalに変更
そしてPick Upオブジェクトを移動すると、グローバルな座標での移動が可能になる
アイテムが1個だと寂しいので、コピーします
Edit > Duplicateからも出来るが、Ctrl＋D（Macは⌘＋D）のショートカットで複製可能
今回は12個配置
アイテムの色変更
Materialを複製して紐付ける
マテリアルの色を黄色に変更してプレハブにセット
プレハブにセットするので、画面内のPickUpの色が全て変更される
- 投稿日：2019-10-20T13:45:25+09:00
C#でフォルダ選択を実装してみた
フォルダ選択を行い、その中のアイテムを一括で取り込むプログラムを実装しました。ワイルドカードを用いることによって特定の拡張子のファイルのみ取り込むこともできます。
構成要素
1.Form1
2.menuStrip
- Menu
-- OpenFolder
3.listBox1
画像１ フォームのデザインフォルダ選択の実装
Form1.csusing System; using System.Windows.Forms; namespace FileSelect { public partial class Form1 : Form { public Form1() { InitializeComponent(); } private void openFolderToolStripMenuItem_Click(object sender, EventArgs e) { FolderBrowserDialog fbd = new FolderBrowserDialog(); //説明文 fbd.Description = "フォルダの選択をしてください。"; fbd.RootFolder = Environment.SpecialFolder.Desktop; //初期値を任意に設定する fbd.SelectedPath = @"C:\Desktop"; //新しいフォルダーの作成ボタン表示の有無 fbd.ShowNewFolderButton = true; if (fbd.ShowDialog(this) == DialogResult.OK) { System.IO.DirectoryInfo di = new System.IO.DirectoryInfo(fbd.SelectedPath); //ワイルドカードでjpgのみ抽出 System.IO.FileInfo[] files = di.GetFiles("*.jpg", System.IO.SearchOption.AllDirectories); foreach (System.IO.FileInfo f in files) { //ファイル名のみlistBox1に表示 listBox1.Items.Add(f.Name); } } } } }実行結果
ファイル選択画面
ファイル名がlistBox1に表示された状態最後に
初心者なので参考程度にお願いいたします。
- 投稿日：2019-10-20T12:54:45+09:00
カメラの視線とメッシュとの交点の座標を数学的に求める
概要
メッシュを張ってカメラでそのメッシュを見たときに、カメラの視線とメッシュとの交点の座標を数学的に計算し表示させるプログラムをUnityで作ります。
結果
さっきの続き
— みこし (@Mikoshi_prog) October 16, 2019
メッシュとカメラの視線が交わるときだけ、交点を表す紫の球体を出現させるようにしました
点の内外判定を使って実装しました?#Unity pic.twitter.com/oa2YpnLkVy方法
まずメッシュをスクリプトを使って張ります。メッシュの張り方はこちらのサイトを参考にさせていただきました。
本記事では、メッシュの頂点に球体を配置し、Unityのシーンからメッシュの頂点の位置を自由に動かせるという方針で実装していくので、参考にさせていただいたスクリプトに多少新たなコードを追加しています。
まず空のゲームオブジェクトに"Mesh_front"という名前を付け(後で"Mesh_back"も作ります)、MeshRenderer、MeshFilter、白いマテリアル、さらに次のスクリプトをアタッチします。
MeshMaker.csusing System.Collections.Generic; using UnityEngine; public class MeshMaker : MonoBehaviour { public GameObject[] sphere; private Vector3[] point; private Mesh mesh; private MeshFilter meshFilter; private void Start() { point = new Vector3[sphere.Length]; for (int i = 0; i < sphere.Length; i++) point[i] = sphere[i].gameObject.transform.position; mesh = new Mesh(); List<Vector3> verticles = new List<Vector3>(); for (int i = 0; i < point.Length; i++) verticles.Add(point[i]); mesh.SetVertices(verticles); List<int> triangles = new List<int>(); for (int i = 0; i < point.Length; i++) triangles.Add(i); mesh.SetTriangles(triangles, 0); meshFilter = GetComponent<MeshFilter>(); meshFilter.mesh = mesh; } private void Update() { } }メッシュの頂点を表す球をsphere、球の頂点の座標をpointとしてあります。ここまでやると、下の図になります。
次に、Mesh_frontに球を渡します。とりあえず三角形のメッシュを作っていくので、球を3つ用意します。3つの球の名前をPoint0, Point1, Point2とし、Scaleを(0.1, 0.1, 0.1)としてあります。Mesh_frontのMakeMesh.csにあるSpheresのSizeを"3"に設定し、Point0, 1, 2を順番に入れます。ここまでやると、下の図になります。
さて、このままだとメッシュは片面にしか張られないので、もう片方の面にも張ります。やり方は簡単です。Mesh_frontを複製し、複製後のものを"Mesh_back"という名前にします。SphereのSizeを"3"にし、球をPoint2, 1, 0の順番に入れます。先ほどとは入れる順番が逆です。
次に、Main Cameraにアタッチするスクリプトについて説明します。FindCoordinateController.csusing UnityEngine; using UnityEngine.UI; public class FindCoordinateOfIntersection : MonoBehaviour { public GameObject[] spheres; //メッシュの頂点の球 public GameObject intersection_sphere; //交点の球 public Text coordinate; //画面上に出す座標のテキスト private Vector3[] point; //メッシュの頂点の座標 private Vector3 normal, intersection; //normal: メッシュの法線ベクトル、intersection: 交点の座標 private MeshMaker meshMaker; private float[] abcd; //メッシュが貼られている平面の方程式の係数 private float parameter; //交点の座標を求めるときに使うパラメータ private bool intersectionIsInsidePolygon; //交点がメッシュの中に存在しているかどうか private void Start() { point = new Vector3[spheres.Length]; for (int i = 0; i < spheres.Length; i++) point[i] = spheres[i].gameObject.transform.position; normal = CalculateOuterProduct(point[0], point[1], point[2]); //3点point[0]~[2]を通る平面の法線ベクトルを求める } private void Update() { abcd = CalculateEquationOfPlane(point[0], point[1], point[2]); //gameObject.transform.rotation * Vector3.forward, gameObject.transform.positionはカメラの視線の方向ベクトル intersection = CalculateCoordinateOfIntersection(abcd, gameObject.transform.rotation * Vector3.forward, gameObject.transform.position); intersectionIsInsidePolygon = WhetherIntersectionIsInsidePolygon(point, intersection, normal); //交点がメッシュの内部にあるときだけ交点をアクティブにする intersection_sphere.SetActive(intersectionIsInsidePolygon); //交点がメッシュの内部にあり、かつ交点がカメラの前側にあるときに、テキストに交点の座標を表示する if (intersectionIsInsidePolygon && parameter > 0f) coordinate.text = "(" + intersection.x.ToString("F2") + ", " + intersection.y.ToString("F2") + ", " + intersection.z.ToString("F2") + ")"; else coordinate.text = "※交点が存在しません※"; } //変数intersectionの値を読み取る、書き込むプロパティ public Vector3 Intersection { get { return intersection; } private set { intersection = value; } } //このメソッドは、vec1,vec2,vec3の3点を通る平面の方程式ax+by+cz+d=0のa,b,c,dを配列で返す private float[] CalculateEquationOfPlane(Vector3 vec1, Vector3 vec2, Vector3 vec3) { float[] ans = new float[]{ normal.x, normal.y, normal.z, -normal.x * vec1.x - normal.y * vec1.y - normal.z * vec1.z }; return ans; } //このメソッドでは、カメラの視線とメッシュとの交点の座標が求められる private Vector3 CalculateCoordinateOfIntersection(float[] plane, Vector3 angle, Vector3 position) { parameter = -(plane[0] * position.x + plane[1] * position.y + plane[2] * position.z + plane[3]) / (plane[0] * angle.x + plane[1] * angle.y + plane[2] * angle.z); float x = angle.x * parameter + position.x; float y = angle.y * parameter + position.y; float z = angle.z * parameter + position.z; return new Vector3(x, y, z); } //このメソッドでは、vec1,vec2,vec3の3点を通る平面の法線ベクトルが求められる private Vector3 CalculateOuterProduct(Vector3 vec1, Vector3 vec2, Vector3 vec3) { Vector3 tmp1 = vec1 - vec2; Vector3 tmp2 = vec1 - vec3; return Vector3.Cross(tmp1, tmp2); //Vector3.Crossは外積を求めるメソッド } //このメソッドは引用させていただきました private bool WhetherIntersectionIsInsidePolygon(Vector3[] vertices, Vector3 intersection, Vector3 normal) { float angle_sum = 0f; for (int i = 0; i < vertices.Length; i++) { Vector3 tmp1 = vertices[i] - intersection; Vector3 tmp2 = vertices[(i + 1) % vertices.Length] - intersection; float angle = Vector3.Angle(tmp1, tmp2); Vector3 cross = Vector3.Cross(tmp1, tmp2); if (Vector3.Dot(cross, normal) < 0) angle *= -1; angle_sum += angle; } angle_sum /= 360f; return Mathf.Abs(angle_sum) >= 0.1f; } }各メソッドの説明をしていきます。
CalculateEquationOfPlaneメソッド
このメソッドの実装はこちらのサイトの1：外積と法線ベクトルを用いる方法を参考にさせていただきました。
p(x-x_0)+q(y-y_0)+r(z-z_0)=0\\ ⇔px+qy+rz-px_0-qy_0-rz_0=0という変形のもとに実装しました。法線ベクトルの値は、すでにStartメソッドの中で計算済みです。
CalculateCoordinateOfIntersectionメソッド
このメソッドの実装はこちらのサイトを参考にさせていただきました。ここで、カメラの座標を表すベクトルを
(p_x, p_y, p_z)カメラの視線の方向ベクトルを
(d_x, d_y, d_z)とします。このときパラメータtを用いて(91)式のように書くとすると
x=d_xt+p_x\\ y=d_yt+p_y\\ z=d_zt+p_zとなります。このx, y, zの値をとる座標がカメラの視線と平面との交点になっていればよいので、平面の方程式を
ax+by+cz+d=0としたとき
a(d_xt+p_x)+b(d_yt+p_y)+c(d_zt+p_z)+d=0\\ ⇔(ad_x+bd_y+cd_z)t=-(ap_x+bp_y+cp_z+d)\\よって左辺が0でなければ
t=-\frac{ap_x+bp_y+cp_z+d}{ad_x+bd_y+cd_z}となります。あとはこのtの値をx, y, zに代入すればOKです。
CalculateOuterProductメソッド
このメソッドはVector3.Crossメソッドを用いて、vec1, vec2, vec3の3点を通る平面の法線ベクトルを求めています。
WhetherIntersectionIsInsidePolygonメソッド
このメソッドはこちらのサイトから、ほぼ完全に引用させていただきました。交点がメッシュの内側にあるのか外側にあるのかを判定します。
メインカメラのInspectorは下の図のようになっています。TPS Cameraというスクリプトは、カメラの位置や回転をXboxコントローラーで制御するためのスクリプトです。本記事では説明を割愛します。また、Find Coordinate Of Intersectionに入っているIntersectionとCoordinateというオブジェクトは後ほど作ります。
次に、カメラとメッシュとの交点を作ります。本記事では交点が視覚的に分かりやすいように、交点の位置に球を配置することにします。球は紫色にし、IntersectionController.csというスクリプトをつけます。
IntersectionController.csusing UnityEngine; public class IntersectionController : MonoBehaviour { private FindCoordinateOfIntersection f; public GameObject camera; private void Start() { f = camera.GetComponent<FindCoordinateOfIntersection>(); } private void Update() { gameObject.transform.position = f.Intersection; } }IntersectionController.csの中では、球の位置を交点の位置に合わせ続けるように実装してあります。Main cameraにFindCoordinateController.csをアタッチしてあるので、Unityでcameraの欄にはMain Cameraを入れておきます。
次に、画面上に交点の座標を表示するためのテキストを作ります。
上の図のようにテキストを配置し、CoordinateというテキストをMain CameraにアタッチしたFindCoordinateController.csに入れておきます。最後に、メッシュが空中に浮いていると頂点の前後関係などが視覚的に分かりづらいと思ったので、x,x,z軸をオブジェクトとして作ることにしました。円柱を細く長く引き伸ばし、軸ごとの色を付けて置いておきます。
以上です。長い記事を読んでくださりありがとうございました！
- 投稿日：2019-10-20T02:44:13+09:00
確率微分方程式の数値解法 (Milstein 法)
Milstein 法
Eular-Maruyama 法は Weak Convergence の次数が 1、Strong Convergence の次数が 1/2 でした。これに Wiener 過程の2次の項を加えることで、Strong Convergence の次数を 1 に上げたものを Milstein 法と呼びます。ノーテーションは前回の記事 と同じです。
Milstein 法の構成
構成法としては、Eular-Maruyama 法では無視した項のうち、Wiener 過程の2次の項を加えます。これに該当する項は下です。
\int_{t_0}^t dW^j_\tau \int_{t_0}^\tau dW_s^k \mathcal{L_k} L^i_j(X_s, s)伊藤の補題を用いて被積分関数を展開し、最低次のみを考慮します。
\int_{t_0}^t dW^j_\tau \int_{t_0}^\tau dW_s^k \mathcal{L_k} L^i_j(X_s, s) \sim \mathcal{L_k} L^i_j(X_{t_0}, t_0) \int_{t_0}^t dW^j_\tau \int_{t_0}^\tau dW_s^kここで問題となるのが Wiener 過程の2重積分です。1次元なら解析的に計算ができますが、多次元の場合は解析的に計算をすることはできません。そのため数値的に積分を行う必要があります。ここを効率的に計算する方法は色々と提案されているようですが (例えば、https://arxiv.org/abs/1107.0151v3 ) ここではもっとも簡単に、以下の方法で数値積分を行いましょう。
\int_{t_0}^t dW^j_\tau \int_{t_0}^\tau dW_s^k = \int_{t_0}^t dW^j_\tau (W_{\tau}^k - W_{t_0}^k) \sim \sum_{l=0}^{N-1} (W_{t_{l+1}}^j - W_{t_l}^j) (W_{t_l}^k - W_{t_0}^k)ここで $t_N = t$ です。
これらを踏まえ、Melstein 法は以下のステップで計算できます。
ステップ1
サンプリング
w^{jk} = \sum_{l=0}^{N-1} (W_{t_{l+1}}^j - W_{t_l}^j) (W_{t_l}^k - W_{t_0}^k)ステップ2
以下の表式を用いて時間発展
\begin{align} t_{i+1} &= t_i +\delta t\\ X_{t_{i+1}}^j &= X_{t_i}^j + f^j(X_{t_i}, t_i)\delta t + L^j_k(X_{t_i}, t_i) w^k + \mathcal{L_k} L^j_l w^{kl} \end{align}まずは各時間グリッドで、各 Wiener 過程のサンプリングを行います。これは例えば以下のようにかけます。
/// <param name="gridStart">t_i</param> /// <param name="gridEnd">t_{i+1}</param> /// <param name="steps">to divide [t_i, t_{i+1}] by this size (decide dt as: dt *steps = t_{i+1} - t_i)</param> /// <returns>matrix[i, j] /// i: index of IEnumerable : time index i [t_0 + dt * i, t_0 + dt * (i+1)] /// j: index of Vector : sample from j'th wiener at the i'th time grid /// </returns> private MathNet.Numerics.LinearAlgebra.Matrix<Complex> ComputeWiernerSamplingMatrix(Parameter gridStart, Parameter gridEnd, int steps) { int dim = SDE.L.Functions[0].Count(); var tmp = Enumerable.Range(0, steps + 1).Select(n => gridStart + n * (gridEnd - gridStart) / steps); var gridEnds = tmp.Skip(1); var gridStarts = tmp.SkipLast(1); var grids = gridStarts.Zip(gridEnds, (start, end) => new { start, end }); var normalByTimeGrid = grids.Select(grid => (new WienerProcess(grid.end, dim)) - new WienerProcess(grid.start, dim)); var randomByTimeGrid = normalByTimeGrid.Select(normal => normal.Generate().Select(r => (Complex)r).ToList()).ToList(); var wienerSamplingMatrix = MathNet.Numerics.LinearAlgebra.Matrix<Complex>.Build.DenseOfRows(randomByTimeGrid); return wienerSamplingMatrix; }このサンプリング結果を用いて Wiener 過程の積分を次のように書くことができます。
/// <summary> /// Compute the \int^t_{i+1}_t_i (W_\tau - W_t_i) dW_\tau /// </summary> /// <param name="intStart">start point of the integral</param> /// <param name="intEnd"> end point of the integral</param> /// <param name="wienerSamplingMatrix">computed by ComputeWiernerRandomValueMatrix</param> /// <param name="integrationSteps">how many steps to integrate the time grid (intEnd - intStart)</param> /// wiener sampling matrix computed by ComputeWiernerRandomValueMatrix /// <returns> /// results matrix[i, j] : \int (W_\tau^i - W_t_0^i) dW_\tau^j</returns> protected Matrix<Complex> ComputeIntegrateWtdWtTerm(Parameter intStart, Parameter intEnd, Matrix<Complex> wienerSamplingMatrix, int integrationSteps) { int wienerDim = SDE.L.Functions[0].Count(); // make time grids for integration var startToEnd = (intEnd - intStart); var dt = startToEnd / integrationSteps; var starts = Enumerable.Range(0, integrationSteps).Select(i => intStart + dt * i); var ends = starts.Select(s => s + dt); var grids = starts.Zip(ends, (start, end) => (start, end)); var wtMinusW0 = wienerSamplingMatrix.ColumnSums(); var wienerIntegralMatrix = Matrix<Complex>.Build.DenseDiagonal(wienerDim, wienerDim, delegate (int ii) { return (wtMinusW0[ii] * wtMinusW0[ii] - startToEnd.Time) / 2; }); for (int i = 0; i < wienerDim; ++i) { for (int j = 0; j < wienerDim; ++j) { if (i == j) continue; Complex tmp = 0.0; for (int k = 0; k < integrationSteps - 1; ++k) { tmp += Enumerable.Range(0, k + 1).Select(l => wienerSamplingMatrix[l, i]).Sum() * (wienerSamplingMatrix[k + 1, j] - wienerSamplingMatrix[k, j]); } wienerIntegralMatrix[i, j] = tmp; } } return wienerIntegralMatrix; }次に、Wiener 過程の積分についている $\mathcal{L_j}L$ を計算します。これは以下のようにかけるでしょう。
/// <summary> /// Compute \mathcal{L}_j L /// </summary> /// <param name="input">input</param> /// <returns> /// return[i][j, k] = \mathcal{L}_i L_jk /// </returns> private IEnumerable<Matrix<Complex>> ComputeMilsteinLjLTerm(Input input) { var dv = MathNet.Numerics.LinearAlgebra.Vector<Complex>.Build.DenseOfEnumerable(((MilsteinSolverCondition)Condition).Infinitesimal); int dim = dv.Count; int wienerDim = SDE.L.Functions[0].Count(); var mathcalLjL = Enumerable.Range(0, dim).Select(i => { var dx = MathNet.Numerics.LinearAlgebra.Vector<Complex>.Build.DenseOfArray(new Complex[dim]); var delta = dv[i]; dx[i] = delta; return ((SDE.L.GetValues(new Input(input.Xt + dx, input.Time)) - SDE.L.GetValues(new Input(input.Xt - dx, input.Time))) / (2 * delta)); }); return mathcalLjL; }上記をまとめることで、Milstein 法で付け加わる項を次のように計算できます。
private MathNet.Numerics.LinearAlgebra.Vector<Complex> ComputeMilsteinTerm(Input input, Parameter intStart, Parameter intEnd, Matrix<Complex> wienerSamplingMatrix) { int dim = input.Xt.Count; int wienerDim = SDE.L.Functions[0].Count(); var differentiateMatrice = ComputeMilsteinLjLTerm(input).ToList(); var lMatrix = SDE.L.GetValues(input); var wIntegral = ComputeIntegrateWtdWtTerm(intStart, intEnd, wienerSamplingMatrix, ((MilsteinSolverCondition)Condition).WienerIntegrationSteps); return MathNet.Numerics.LinearAlgebra.Vector<Complex>.Build.DenseOfEnumerable(Enumerable.Range(0, dim).Select(i => { var matrix = Matrix<Complex>.Build.Dense(dim, wienerDim); // to compute with Trace, arrange the index for (int j = 0; j < dim; ++j) { for (int k = 0; k < wienerDim; ++k) { matrix[j, k] = differentiateMatrice[j][i, k]; } } return (matrix * wIntegral * lMatrix.Transpose()).Trace(); })); }上で計算されたもの以外は Eular-Maruyama 法と同様に計算できるため、Eualar-Maruyama 法のクラスを継承し、ランダム項を計算する関数をオーバーライドすれば十分です。
protected override MathNet.Numerics.LinearAlgebra.Vector<Complex> ComputeRandomTerm(Input input, Parameter start, Parameter end) { var randomMatrix = ComputeWiernerSamplingMatrix(start, end, ((MilsteinSolverCondition)Condition).WienerIntegrationSteps); var lMatrix = SDE.L.GetValues(input); var randomVector = randomMatrix.ColumnSums(); var eularMaruyamaTerm = lMatrix * randomVector; var milsteinTerm = ComputeMilsteinTerm(input, start, end, randomMatrix); var randomTerm = eularMaruyamaTerm + milsteinTerm; return randomTerm; }実験
せっかくなので多次元確率微分方程式をときましょう。ここでは以下の微分方程式を解くことを考えます。
\begin{align} dX_t &= rX_t dt + \sigma_1 X_t dW^1_t\\ dr_t &= a(b-r_t)dt + \sigma_2 dW^2_t + \sigma_2 dW^1_t \end{align}1つ目の方程式は BS モデルであり、2つ目は Vasicek モデル (https://ja.wikipedia.org/wiki/バシチェック・モデル) に Wiener 過程の項を1つ付け加えたものです。Vasicek モデルに Wiener 項をひとつ付け加えた理由としては、積分計算の際に作る行列の非対角成分がゼロとならないようにするためです。
Vasicek モデルは金融における短期金利のモデルの1つであり、Hull-White モデル (https://ja.wikipedia.org/wiki/ハル・ホワイト・モデル) の一部パラメータを定数としたものです。Hull-White モデルは、BS モデルのように対数正規分布ではなく、マイナス金利を表現できるのが強みです (詳しくは Shreve や Bjork などの数理ファイナンス本)。Vasicek モデルは解析的に解くことができ、今回の2つ目の方程式に書き換えても解析的に解くことができます。
\begin{align} r_t &= r_{t_0} e^{-a(t-t_0)} + b(1-e^{-a(t-t_0)}) + \sigma_2 e^{-a(t-t_0)} \int_{t_0}^t e^{as} dW^1_s + \sigma_2 e^{-a(t-t_0)} \int_{t_0}^t e^{as} dW^2_s\\ \mbox{E}[r_t] &= r_{t_0} e^{-a(t-t_0)} + b(1-e^{-a(t-t_0)})\\ \mbox{V}[r_t] &= \frac{\sigma_2^2}{a} \end{align}ここでは以下の条件の下で実験を行いました。
\begin{align} t_0 &= 0\\ t &= 1\\\\ r &= 0.01\\ \sigma_1 &= 0.1\\ X_0 &= 1\\\\ a &= 0.1\\ b &= 0.1\\ \sigma_2 &= 0.01\\ r_0 &= 0 \end{align}確率微分方程式を解く際の時間グリッドは等間隔で100個に分け、Wiener 過程の積分に関しても等間隔に、10 のグリッドとしました。$\mathcal{L_j}$ の微分を計算する際は数値的に微分を行い、微分に用いる微小量は $10^{-10}$ とし、中央差分で数値微分を行いました。
発生させたパス数は 100000 で、それらパスから標本平均、標本分散を計算しました。数値計算の結果は以下となりました。良い精度で計算ができていることがわかります。
$X_t$ の数値計算結果
Milstein 法 解析解 平均 $1.0099$ $1.0101$ 分散 $9.3879 \times 10^{-3}$ $1.0253 \times 10^{-2}$ $r_t$ の数値計算結果
Milstein 法 解析解 平均 $9.5128 \times 10^{-3}$ $9.5163 \times 10^{-3}$ 分散 $1.7911 \times 10^{-4}$ $1.8127 \times 10^{-4}$