雑感等

音楽,数学,語学,その他に関するメモを記す.

文献の示し方

出典:若林敦. 理工系の日本語作文トレーニング. 朝倉書店, 2000.

(i) 図書を示す場合([ ] は書かない場合もある):

著(編)者名 書名 [版表示] 出版社(者)またはシリーズ記載 [出版地] 出版年 [ページ]

[版表示]=初版は不要.改訂版・増補版及び第2版以後の版は必ず表示する.
[ページ]=引用箇所,参照箇所を特に指定する場合に表示する.
[出版地]=日本語の文章では,書かないことが多い.

(ii) 図書収蔵論文を示す場合:

論文の著者名 論文名(表題) 所収図書の編(著)者名 書名 …途中は,上の図書を示す場合に同じ… [ページ]

[ページ]=記載する場合は所収論文の掲載されているはじめのページとおわりのページを記す.おわりのページは省略してもよい.

(iii) 雑誌論文,雑誌記事を示す場合:

著者名 記事タイトル,論文名(表題) 誌名 巻(号) 発光年 [ページ]

[ページ]=記載する場合は論文,記事の掲載されているはじめのページとおわりのページを記す.おわりのページは省略してもよい.

フラクショナルブラウン運動の生成するmatlabプログラム(Vossのアルゴリズムによる)

function sig=voss(hurst,len)
%Vossのアルゴリズムによってフラクショナルブラウン運動を生成する.
%husrt:Hurst指数
%len:出力するデータ長

%2^n+1>=lenとなる最小のnを求め,その長さでfBmを生成する
lenFbm=2^ceil(log(len-1)/log(2))+1;

%初期化
fBm=zeros(lenFbm,1);

%fBmの終点と始点の値を決める.
fBm(lenFbm)=0;
s2_0=1; %fBm(lenFbm)の値を決める乱数の分散.このパラメタによりfBm全体の振幅が変わる.
fBm(lenFbm)=randn()*sqrt(s2_0);

%再帰的に残りの点を決める.
fBm=calcmidpoint(1,1,lenFbm,fBm,hurst,s2_0);

%出力するデータの長さを指定された長さで打ち切る
sig=fBm(1:len);

end

function fBm=calcmidpoint(n,a,b,fBm,hurst,s2_0)
%n:再帰の深さ
%a:区間の左端要素の添字
%b:区間の右端要素の添字
%fBm:フラクショナルブラウン運動の時系列データ
%hurst:ハースト指数
%s2_0:fBm(end)の値を決めるのに使った乱数の分散

%中点の添字を計算
mid=(b-a)/2+a;

%中点が整数番目の添字を持つなら中点の値を計算する.
if mid==fix(mid)
    
    %s2:分散
    s2=( ( 1-2^(2*hurst-2) )*s2_0 )/( 2^(2*hurst*n) );
    
    %中点の値=区間の両端の値の平均値+分散
    fBm(mid)=(fBm(b)+fBm(a))/2+randn()*sqrt(s2);
    fBm=calcmidpoint(n+1,a,mid,fBm,hurst,s2_0);
    fBm=calcmidpoint(n+1,mid,b,fBm,hurst,s2_0);
end
end

Vossのアルゴリズムに関する式変形の過程(フラクショナルブラウン運動の生成)

文献: 本田勝也. フラクタル. 朝倉出版, 2002. pp. 110-112.

ラクショナルブラウン運動の関数を X(x),\, (0\leq x \leq1)を生成する「ヴォスのアルゴリズム」が文献に示されていた.
その中の式変形を詳細に示す.

 X\left(  0 \right) =0,\, X\left(  1 \right) =\delta_{0}, \, (\delta_{0} \sim \text{平均0,分散1の正規分布})
 X \left( \frac{1}{2} \right) = \frac{1}{2} \left\{ X(0)+X(1) \right\} + \delta_{1},\,\left( \delta_{1} \sim \text{平均0,分散$\sigma_{1}^{2}$の正規分布} \right) \hspace{3em}\text{(10.17)}
ただし X(1)-X(0) \delta_1は互いに独立である.

 \langle \left\{  X \left( \frac{1}{2} \right) - X \left(  0 \right)   \right\}^{2} \rangle \hspace{3em}\text{(10.18)左辺}
式(10.17)を代入
 =\langle \left\{  \frac{1}{2} \left\{ X(0)+X(1) \right\} + \delta_{1} - X \left(  0 \right)   \right\}^{2} \rangle
展開
 =\langle \left\{  \frac{1}{2}  X(0) + \frac{1}{2} X(1)  + \delta_{1} - X \left(  0 \right)   \right\}^{2} \rangle
項をまとめる
 =\langle \left\{   \frac{1}{2} X(1) -\frac{1}{2}  X(0)  + \delta_{1}    \right\}^{2} \rangle
 =\langle \left\{   \frac{1}{2} \left\{ X(1) - X(0) \right\} + \delta_{1}    \right\}^{2} \rangle
二乗を展開
 =\langle   \frac{1}{4} \left\{ X(1) - X(0) \right\}^{2} + {\delta_{1}}^{2} +  {\delta_{1}} \left\{ X(1) - X(0) \right\} \rangle
加法性を適用(統計平均と統計の期待値(平均)は同じ?)
 =\langle  \frac{1}{4} \left\{ X(1) - X(0) \right\}^{2} \rangle + \langle {\delta_{1}}^{2} \rangle + \langle  {\delta_{1}} \left\{ X(1) - X(0) \right\} \rangle
 X(1)-X(0) \delta_1は互いに独立だから,積の平均を平均の積に分解
 =\langle  \frac{1}{4} \left\{ X(1) - X(0) \right\}^{2} \rangle + \langle {\delta_{1}}^{2} \rangle + \langle  {\delta_{1}}\rangle  \langle  \ X(1) - X(0) \rangle
正規分布に従う乱数の平均だから   \langle {\delta_{1}}^{2} \rangle = \sigma_{1}^{2}, \,   \langle  {\delta_{1}} \rangle= 0
 =\langle  \frac{1}{4} \left\{ X(1) - X(0) \right\}^{2} \rangle + \sigma_{1}^{2} + 0 \times  \langle  X(1) - X(0) \rangle
 =\langle  \frac{1}{4} \left\{ X(1) - X(0) \right\}^{2} \rangle + \sigma_{1}^{2} \hspace{3em}\text{(10.18)右辺}


他の書籍を読んだ際,統計平均を\langle  \rangleで表したり,ほぼ等しいことを \sim で表したりする記法を知らなかったため,式が理解できないことがあった.
専門分野の表記はよそ者がみると理解できない.

フラクショナルブラウン運動の分散特性のプロット

matlabプログラムで分散特性\sigma^{2}_{H}= \langle \lvert f(t+\tau)-f(t) \rvert ^{2} \rangleに対する時間間隔 \tauを両対数グラフにプロットする.

Hurst指数がHのフラクショナルブラウン運動(fBm)の信号に対して分散 \sigma^{2}_{H}-時間間隔 \tauをプロットすると,
両対数グラフ上で傾きが2Hの線形となる範囲が現れる.

フラクショナルブラウン運動の生成プログラム - 雑感等に記載したプログラムを用い,
Hurst指数H=0.2, 0.5, 0.8のfBm信号を生成し,分散 \sigma^{2}_{H}-時間間隔 \tauをプロットした.

 \tauが小さい領域では線形にプロットされる.

H=0.2
f:id:kazmus:20180804111504p:plain

H=0.5
f:id:kazmus:20180804111517p:plain

H=0.8
f:id:kazmus:20180804111529p:plain

分散 \sigma^{2}_{H}-時間間隔 \tauをプロットするmatlabプログラム

%使用例
%fBm=csvread('fBm.csv');
%dif_time(fBm)
%第二引数は空で良い

function dif_time(sig,fs)
if 1~=exist('fs')
    fs=length(sig);
end
ts=1/fs;
loopMax=length(sig);
meanMoment=zeros(loopMax,1);
subplot(211)
for k=1:100:loopMax
    diffSquareSig=(sig(k+1:end)-sig(1:end-k)).^2;
    loglog((k*ts), ( mean(diffSquareSig)) ,'o')
    hold on
end
ylabel('<|f (t+\tau)-f (t)|^{2}>','FontSize',16)
xlabel('\tau','FontSize',16)
subplot(212)
plot(sig)
ylabel('amplitude','FontSize',16)
xlabel('time','FontSize',16)
set(gca,'FontSize',12)
end

フラクショナルブラウン運動の生成プログラム

正規分布に従う乱数を分数階積分してフラクショナルブラウン運動を生成し,csv形式に出力するプログラム.
下記のC++11のコンパイラで動作を確認

  • Embarcadero Technologies Inc. bcc32c version 3.3.1
  • Microsoft Visual Studio Community 2015 Version 14.0.25431.01 Update 3
//>>bcc32c fBm_generator.cpp
//>>fBm_generator.exe 
//または
//>>fBm_generator.exe [Hurst指数]
#include<iostream>
#include<math.h>
#include <random>
#include <fstream>
#include <stdlib.h>

////////設定項目////////

//Hurst指数 下記の値または,コマンドライン引数の一番目
#define HURST_EXPONENT (0.5)
//出力する時系列データの点数
#define DATALEN (65536)

//fBmの時系列データを出力するファイル名.ただし,出力形式はcsv
#define fBm_FILENAME ("fBm.csv")

//fBmの生成に使用する乱数の時系列データを出力するファイル名.ただし,出力形式はcsv
#define WHITE_NOISE_FILENAME ("white.csv")


////////計算プログラム////////

//分数階積分する関数"FractionalIntegral"は下記のサイトから引用,一部改変した.
//http://www.geocities.co.jp/Technopolis-Mars/1795/study/fractint.html
#define MAXKAISUU s //積分演算をおこなう時に利用する過去の情報の最大数
#define m_chk_mean false //最終的に平均値を零とする処理をおこなうかどうか
#define m_DT 1.0 //時系列のサンプリング時間間隔

//入力信号をそのまま微分する
//*data・・・微分をおこなう対象データのポインタ
//sekibun・・・積分階数
//s・・・データの大きさ
void FractionalIntegral(double *data, double sekibun, int s)
{
	register int m, t;

	//sekibunが整数かどうか判別する(整数の時には動作が異なる)
	int kaisuu = (sekibun == (int)sekibun) ? (int)sekibun : MAXKAISUU;
	double* ganma = new double[kaisuu + 1];

	//ガンマ関数の値をストックする    (n+1)C(m+1)
	//ついでに符号も保存する
	ganma[0] = 1.0;
	for (m = 1; m <= kaisuu; m++)
		ganma[m] = -ganma[m - 1] * (sekibun - m + 1) / (double)m;

	//非整数階積分をする
	double p = pow(m_DT, sekibun); //h^d
								   // for(t = s; t; data[--t] -= data[0]);  //積分の場合はいらない

	register double f;
	int M;
	for (t = 0; t < s; t++) {
		M = (t < kaisuu) ? t : kaisuu;
		for (f = data[t] * p, m = 1; m <= M; f -= ganma[m] * data[t - m], m++);
		data[t] = f;
	}
	delete[] ganma;

	//データの平均を0に
	if (m_chk_mean) {
		for (f = 0.0, t = s; t; f += data[--t]);
		for (f /= s, t = s; t; data[--t] -= f);
	}
}

//正規分布に従う乱数を発生させる関数"randnormal"は下記のサイトを引用,改変した.
//https://cpprefjp.github.io/reference/random.html
double randnormal(double mean, double stdDev) {
	// メルセンヌ・ツイスター法による擬似乱数生成器を、
	// ハードウェア乱数をシードにして初期化
	std::random_device seed_gen;
	std::mt19937 engine(seed_gen());

	// 正規分布
	// 平均mean、標準偏差stdDevで分布させる
	std::normal_distribution<> dist(mean, stdDev);
	return dist(engine);
}


int main(int argc, char *argv[]) {
	double signal[DATALEN];

	//正規分布に従う乱数を発生
	for (int i = 0; i < DATALEN; ++i) {
		signal[i] = randnormal(0, 1);
	}
	std::cout << "乱数の一番目の値:" << signal[0] << std::endl;
	//正規分布に従う乱数を保存
	std::ofstream filew(WHITE_NOISE_FILENAME);
	for (int i = 0; i < DATALEN; ++i) {
		filew << signal[i] << "\n";
	}

	//Hurst指数に0.5を加算し,乱数時系列信号を分数階積分する階数integrateOrderを決定
	double integrateOrder = HURST_EXPONENT + 0.5;

	if (argc > 1) {
		char *strtod_error;
		double input = strtod(argv[1], &strtod_error);
		if (*strtod_error != '\0') {
			std::cout << "Hurst指数として設定できない文字列:" << strtod_error << std::endl;
		}
		integrateOrder = input + 0.5;
	}
	std::cout << "Hurst指数:" << integrateOrder - 0.5 << std::endl;
	//分数階積分を実行
	FractionalIntegral(signal, integrateOrder, DATALEN);

	//フラクショナルブラウン運動信号を保存
	std::ofstream file(fBm_FILENAME);
	for (int i = 0; i < DATALEN; ++i) {
		file << signal[i] << "\n";
	}

	file.close();
	filew.close();
	return 0;
}

指数関数表記のシフト演算子

(p. 174)方程式 x' (\lambda)=wx(\lambda)を満たし,かつ条件 x(\lambda_{0})=kを満たすような演算子関数 x(\lambda)はたかだか一つしか存在しない.
(p. 175) wが数ならば,指数関数 x(\lambda)=e^{\lambda w}は方程式 x' (\lambda)=wx(\lambda)および条件 x(0)=1を満足する.
(p. 172) h^{\lambda}は, (h^{\lambda})'=-sh^{\lambda},\, h^{0}=1を満たす.
(p. 175)シフト演算子(移動演算子 h^{\lambda}=e^{-\lambda s}と書く.

  1. ヤン・ミクシンスキー. 演算子法: 上巻. 松村英之, 松浦重武訳. 新版, 裳華房, 299p.

不変測度(不変密度)がフロベニウス=ペロン方程式を満たすことの証明

文献[1]でロジスティック写像リアプノフ指数を解析的に求めるために,不変測度とフロベニウス=ペロン方程式について述べていた.

  • ロジスティック写像 x_{n+1}=f(x_n)=4 x_{n} (1-x_{n})

以下では不変測度とフロベニウス=ペロン方程式の詳細についてまとめる.

不変測度(不変密度)がフロベニウス=ペロン方程式を満たすことの証明

  • 不変測度: \rho(x)=\lim_{T \rightarrow \infty} {\frac{1}{T} \sum^{T-1}_{n=0} {\delta (x-x_{n})}}
  • フロベニウス=ペロン方程式: \rho(x)= \int d\xi \delta (x - f(\xi)) \rho (\xi)
  • フロベニウス=ペロン方程式に不変測度を代入する.

 \rho(x)= \int d\xi \delta (x - f(\xi)) \{ \lim_{T \rightarrow \infty} {\frac{1}{T} \sum^{T-1}_{n=0} {\delta (\xi-x_{n})}} \}

  \int d\xi \delta (\xi - f(\xi))は極限の Tに依存しないから,  \int d\xi \delta (x - f(\xi))を移動する.
 \rho(x)= \lim_{T \rightarrow \infty} \frac{1}{T} \int d\xi \delta (x - f(\xi))   \sum^{T-1}_{n=0} {\delta (\xi-x_{n})}

  • 積分と総和の順番を変える.

 \rho(x)= \lim_{T \rightarrow \infty} \frac{1}{T}  \sum^{T-1}_{n=0} \int d\xi \delta (x - f(\xi))   \delta (\xi-x_{n})

 \rho(x)= \lim_{T \rightarrow \infty} \frac{1}{T}  \sum^{T-1}_{n=0} \delta (x - f(x_{n}))

  •  x_{n+1}=f(x_n)を用いる.

 \rho(x)= \lim_{T \rightarrow \infty} \frac{1}{T}  \sum^{T-1}_{n=0} \delta (x - x_{n+1})

上式までの変形から, x_nで表される不変測度 \rho (x)をフロベニウス=ペロン方程式に代入すると,不変測度 \rho (x) x_{n+1}で表されることが示される.

 x_n=x_{n+1}が成立しなければ不変測度 \rho (x)をフロベニウス=ペロン方程式を満たさない.不変測度 \rho (x)をフロベニウス=ペロン方程式を満たすためには, x_nが「定常状態」であれば良い.

  • 従って x_n=x_{n+1}が成立すると仮定する.

 \rho(x)= \lim_{T \rightarrow \infty} \frac{1}{T}  \sum^{T-1}_{n=0} \delta (x - x_{n})

以上から,不変測度(不変密度)がフロベニウス=ペロン方程式に対して不変であることを示した.

不変測度(不変密度)

文献[2]では,「軌道上の状態 xが微小区間 [  x, x+dx ]に滞在する割合が \rho (x) dxと表させるとき, \rho (x)を不変密度(invariant density)という.」と述べている.

文献[1], [2]ではエルゴード仮説が正しいと仮定して,時間平均と統計平均(アンサンブル平均)とが等しいと仮定している.

  • 関数 h (x)の時間平均: \lim_{T \rightarrow \infty} {\frac{1}{T} \sum^{T-1}_{n=0} {h( x_n )} }
  • 関数 h (x)の統計平均:  < h > = \int {\rho (x) h (x) dx}

フロベニウス=ペロン方程式

文献[2]では,おおよそ次のように説明している.

初期時刻において滑らかな確率密度関数 \rho_0従う無数の初期点を写像 x_{n+1}=f(x)で発展させる.このとき時刻 nにおいて微小区間 [  x, x+dx ]に存在する点の割合が \rho_n (x) dxと表されるとする.この \rho_n (x)の時間発展は以下のフロベニウス=ペロン方程式に従う.

  • フロベニウス=ペロン方程式: \rho_{n}=\int \rho_n (y) \delta \{ x-f(y) \} dy

さらに同書には,「ここで右辺は,時刻 n確率密度関数 \rho_n (y)に従って分布する状態点 yのうち
写像 fにより新たな状態点 x= f(y)に移される点のみを集めたものが,時刻 n+1における状態 xの密度 \rho_{n+1} (y)を与えることを意味している.」と示されている.

文献

[1] 中川匡弘: カオス・フラクタル感性情報工学. 日刊工業新聞社, 2010.
[2] 中尾裕也ほか: ネットワーク・カオス―非線形ダイナミクス複雑系と情報ネットワーク. コロナ社, 2018.