Symbolic mathで微分演算子っぽいものを作る

Symbolic mathで微分演算子っぽいものを作ることができる。なんの役に立つか不明だが、面白い機能なので覚書。

 

[問題の定式化]

2回微分可能な関数 u(x)について以下の簡単な常微分方程式を考える。

 \frac{d^2 u }{d x^2} + \rho u + x = 0

微分演算子 Dを以下のように定義する。

 D =  \frac{d^2 }{d x^2} + \rho

微分演算子 Dを使えば、常微分方程式は以下のように簡単に表記できる。

 D u + x = 0 

 

関数の定義:

function [DD] = diffOperator(funcArray, vars, params)

rho = params.rho;

xx = vars.xx;

DD = diff(funcArray, xx, 2) + rho * funcArray;

end

 

使い方:

% 関数を特定した場合の使い方

syms xx;

params.rho = 1.5;

vars.xx = xx;

uu(xx) = xx^2; 

func = symfun(uu(xx), xx);

diffOperator(func, vars, params)

 

% 関数が一般形の場合の使い方

clear uu;

syms uu(xx);

func = symfun(uu(xx), xx);

 

diffOperator(func, vars, params)

 

 

Symbolic mathではシンボリック関数の配列が作成できる。微分もできる。

使い方次第だと思うが、面白い機能なので覚書。

 

シンボリック関数はsymfunで作ることができる(フォーマルなやり方)。簡易なやり方としてsym f(x,y)のようにすることもできる(xxとyyは自動生成される)。

 

syms phi1(xx) phi2(xx);

phi1(xx) = xx;
phi2(xx) = xx^2;

funcArray = symfun([1, phi1(xx), phi2(xx)], xx);

 

funcArray(xx) = [ 1, xx, xx^2] 

 

使い方:funcArray(2)とするとans = [ 1, 2, 4]が返ってくる。

 

補足:2変数関数にもできる

clear;

syms phi1(xx, yy) phi2(xx, yy);  % xxとyyは自動生成される

phi1(xx, yy) = xx*yy;
phi2(xx, yy) = xx^2 * yy^2;

funcArray = symfun([1, phi1(xx, yy), phi2(xx, yy)], [xx, yy]);

 

funcArray(xx, yy) = [ 1, xx*yy, xx^2*yy^2] 

 

 使い方:funcArray(2, 2)とするとans = [ 1, 4, 16]が返ってくる。

 

さらにその「シンボリック関数の配列」を(まるで普通の関数のように)微分することもできる。

chebVec(x, y) = chebyshevT([0, 1, 2], x) .* chebyshevT([0, 1, 2], y)

chebVec(x, y)  = [ 1, x*y, (2*x^2 - 1)*(2*y^2 - 1)]

diff(chebVec, x)

ans(x, y) = [ 0, y, 4*x*(2*y^2 - 1)]

 

公式ページのマニュアルだと少しわかりにくい。

参考:

jp.mathworks.com

jp.mathworks.com

 

シンボリック式・シンボリック関数はstringで文字列に変換できる

何に役立つのか分からないが、面白い機能なので覚書(この機能はMatlabドキュメントに載っていないようなので)。

 

syms xx yy;

% シンボリック式はstringで文字列に変換できる
string(xx + yy)

% シンボリック関数もstringで文字列に変換できる
testFun = @(xx, yy) xx + yy;

string(testFun(xx, yy))

 

なお、普通のMatlab関数をstringで文字列に変換・・・はできない模様。

Matlabのシンボリック式を無名関数(関数ハンドル)に変換する関数

Matlabのsymbolic math toolboxを使わない人にとってはどうでもいい話題ですが、シンボリック式を最適化しようとして少し困ったので、覚書。

 

Matlabのシンボリック式またはシンボリック関数を無名関数(関数ハンドル)に変換する関数が存在する。使い方はこんな感じ(↓)

ht = matlabFunction(r)

 

詳しくは

https://jp.mathworks.com/help/symbolic/matlabfunction.html

Miranda-Fackler (2002) p. 132のテンソル積による計算結果がMatlabのkronの結果と同一であることの確認。

% Miranda-Fackler (2002) p. 132のテンル積による計算結果がMatlabのkronの結果と同一であることの確認。

 

syms x1 x2;

aa = [ 1, x1, x1^2]

bb = [ 1, x2]

kron(bb, aa)

 

注:ベクトルの場合、テンソル積とクロネッカー積の結果は同一となる(だったはず)。テンソルは20年位やっておらず、知識が極めて怪しいので、正確に知りたい人はテンソル解析の教科書を見てください。

 

[参考文献]

reference.wolfram.com

MacOSでAnaconda (python)のバージョンアップ

condaを使う場合、

$ conda update conda
$ conda update anaconda
$ conda update --all

Anaconda Navigatorを起動して、バージョンが最新版になっていることを確認する。

TensorFlowとKerasのインストール

$ conda install -c jjhelmus tensorflow
$ conda install -c conda-forge keras

RでX-13ARIMA-SEATSを使う(2019年6月版)

自分用の備忘録です。CRANにあるx13binaryが使えるようになったので、自分でX-13ARIMA-SEATSをコンパイルする必要がなくなりました。

[インストール]GNU Rで二つのパッケージをインストール

  1. install.packages("seasonal"); install.packages("x13binary")

[GNU Rのパッケージseasonalからx13asを実行する]

  1. パッケージの読み込み
    library(x13binary); library(seasonal)
  2. checkX13()を実行して、"Congratulations! 'seasonal' should work fine!"と表示されれば問題なし
  3. m <- seas(AirPassengers)
  4. plot(m)
  5. 季節調整済み系列(s11)の取り出し
    m$series$s11
  6. seasonalの使い方についてvignettesを参照されたい
    https://cran.ism.ac.jp/web/packages/seasonal/vignettes/seas.pdf

[注意すべき点]

  1. X-13ARIMA-SEATSから季節調整済み系列はs11となった模様(X-12-ARIMAのd11。プログラムする人はお気をつけ下さい)
    http://www.census.gov/ts/x13as/docX13ASHTML.pdf

[s10やs11等のまとめ]
s10: 季節成分(seasonal component)
s11: 季節調整値(seasonal adjustment component)
s12: トレンド成分(trend component)
s13: 不規則成分(rregular component)
詳しくは以下のTable7.13を参照:
https://www.census.gov/ts/x13as/docX13AS.pdf