HOME > Ox

Ox

よく使う関数

※詳細はヘルプを参照してください.OxEditの場合,helpから参照可能です.

一般的な関数

rows 行列の行数.See columns.
columns 行列の列数.See rows.
sizeof 配列の次元を取得.2次元配列の場合,行数×列数.

行列関係

diag 対角行列の作成.
diagcat 2つの行列の対角要素で結合.非対角要素はゼロになる.
diagonal 行列の対角要素をベクトル化.
diagonalize 行列の非対角要素をゼロにして,対角行列へ変更.
dropc, dropr 指定したindexの列,行を削除.
deletec, deleter 指定した値を含む列,行を削除.
unit(n) n×nの単位行列の作成.
ones(m,n) m×nの要素が全て1の行列を作成.
zeros(m,n) m×nの要素が全て0の行列を作成.

入力・出力

loadmat 行列データの読み込み.
savemat 行列データの保存.
fopen ファイルを開く.
fclose ファイルを閉じる.
fprint ファイルに書き込む.

繰り返し文

do{ ... } while() 条件を満たす間は繰り返し計算を続ける.条件は計算後に判定.
for 指定回数のループ.
while 条件を満たす間は繰り返し計算を続ける.計算をする前に条件判定を行う.

条件文

if()
switch()

その他

array 配列.
double 倍精度に型変換.
fabs 絶対値を返す.
format 書式を設定.(1)数値の表示形式と(2)一行辺りの表示文字数の変更が可能.

(Uploaded: 2010年10月31日, Last Updated: 2011年03月14日)  ページ上部へ移動

Oxで計量経済学(初級編)

Oxで計量分析をする時のポイントを整理しています.パッケージのように全てを自動で処理してくれるわけではありませんので,自分で処理する必要があります.

変数の宣言

decl x;

使用する変数は必ず事前に宣言declをしなければなりません.上記では変数xを宣言しています.もう一つ注意することは,コマンドの最後にはセミコロンをつけて閉じることです.

データの読み込み

decl data = loadmat("data.csv");

数値データを行列として読み込みとき,loadmat()を使います.上記では,data.csvファイルを読み込み,変数dataに格納しています.データファイル上の文字列は読み込まれません.第一行は特別に変数名として読み込むことが可能です.読込可能な対応ファイルとしてよく使っているのは,csvとxls(エクセル形式)の2つです.

欠損値を除く

deleter( data );

Stata等のパッケージでは,欠損値は自動的に除去されて計算処理されます.Oxでは手動で欠損値を削除する必要があります.deleter()はNaNを含む行を削除します.

特定の値を持つ行を除く

deleter( data, 1 );
deleteifr( data, data[][1] .== 1 );

ある値を持つ行を削除したい場合にdeleter()もしくはdeleteifr()を使います.Stataではdrop ... if ... に対応します.

最小二乗法(OLS)

decl vb, mIXX, mXX;
olsc( vy, mx, &vb, &mIXX, &mXX );

OLS用の関数olsc()が用意されています.olscは得られた推定値を列ベクトルとして返します.また,(X'X)-1X'Xを返すこともできます.上記ではそれぞれ,mIXX, mXXとして変数を宣言しています.vbがOLS推定値です.あとは,アドレス演算子&を使い,事前に宣言した変数にolscからの計算結果を格納します.

結果の出力

結果を表示するにはprint()もしくはprintln()を使用します.println()は文の最後に改行が入ります.

println("vb = ", vb);

表示を整えられるようになっています.例えば,以下では2×2行列mxの列と行にラベルをつけています.その他の詳細はヘルプを参照してください.

println("%c", {"c1", "c2"},
     "%r", {"r1", "r2"}, mx);

ベクトルの内積でスカラーを得たとき,型が行列になってしまうので,スカラーへ変換します.例えば,倍精度に変換する場合,double()を使います.

double( vx'vx );

(Uploaded: 2011年04月15日, Last Updated: 2011年06月09日)  ページ上部へ移動

Oxで計量経済学(中級編)

関数の作成

fMyfun( ){

}

オリジナルの関数を作成し呼び出すこともできます.

パッケージ化

作成中
class MyOLS{
  
  MyOLS();
  LoadVarNames();
  LoadData();
  DoEstimation();
  Estimate();
  Covar();
  CalcStat();
  Output();
  Predict();
  
};
MyOLS::MyOLS(){

}
MyOLS::LoadVarNames(){

}
MyOLS::LoadData(){

}
MyOLS::DoEstimation(){

}
MyOLS::Estimate(){

}
MyOLS::Covar(){

}
MyOLS::CalcStat(){

}
MyOLS::Output(){

}
MyOLS::Predict(){

}

よく使うものはクラスとしてまとめておくと便利です.グローバル変数を使用しなくて済みます.

(Uploaded: 2011年04月15日, Last Updated: 2011年06月17日)  ページ上部へ移動

いろいろなTips

連番ファイル名

データ加工の時に連番のファイル名をつけたい時があります.数値を文字列に変換するためにsprint関数を使用します.

decl sFile;
decl mX;
mX = <0.8, 0.2; 0.1, 0.9>
for(decl i=1; i<10; ++i){
  sFile = "Matrix" + sprint(i+1) + "csv"
  mX = mX^(i+1);
  savemat(sFile, mX);
}

NaN

時系列でラグや階差を取りたい時があります.そのようなときlag0diff0を使用しますが,欠損値はデフォルトではゼロになってしまいます.その代わりにNaNを代入することで,数値データとは区別します.以下はlag0関数の例です.

#include <oxfloat.h> //M_NANのために必要

decl vx, vxl1, vxl2;
vx = rann(10,1);
vxl1 = lag0( vx, 1, M_NAN );
vxl2 = lag0( vx, 2, M_NAN );

もし,上記の3つのベクトルをデータとして使用したいときは,deleter関数でNaNを含む行を削除します.

decl mx;
mx = deleter( vx ˜ vxl1 ˜ vxl2 );

図でギリシャ文字出力

oxdraw.hを読み込むことでグラフへの出力が可能になります.文字列としてLaTeXスタイルの文字列を使えます.ギリシャ文字は以下のようにタイプします.例はβです.

$¥¥beta$

クロネカー積

クロネカー積は**演算子で計算できます.

decl mA, mB, mK;
mK = mA ** mB;

ダミー変数の作成

コードに従ってダミーを作成します.

decl cN;
decl mDummy;
decl temp;
cN = rows(data);
mDummy = zeros(cN,47);

for(decl i=0; i<47; ++i){
  temp = isdotfeq(data[][1],i+1);
  mDummy[][i] = temp;
}

大きな整数の掛け算

計算結果がおかしい??

decl dx = 1400;
println( "%15.5g", 1/(dx^3) ~ 1/(dx*dx^2) ~ 1/(dx*dx*dx) ~ (1/dx)*(1/dx)*(1/dx) );

3.6443e-010  -6.4476e-010  -6.4476e-010  3.6443e-010  3.6443e-010

もしこのまま計算し続けると,符号自体も変わってしまうので注意.なお,これと同じ計算をRやExcelでやっても標準ではこのようなことは起こりません.

理由.int型の上限値を超えているから.大きい正の整数の掛け算をすると負になる.大きな整数で階上のような計算をそのまま打ち込むと間違えます.

解決策.int型で保存されているのが問題なので,double型で保存すれば解決できます.整数の計算が大きくなる時は一つでもint型をdouble型をしていおけば,計算結果はdouble型で返されます.例えば,1400.に書き換える等すれば可能です(1400の後ろに点が付いていることに注意).

decl dx = 1400.;
println( "%15.5g", 1/(dx^3) ~ 1/(dx*dx^2) ~ 1/(dx*dx*dx) ~ (1/dx)*(1/dx)*(1/dx) );

3.6443e-010  3.6443e-010  3.6443e-010  3.6443e-010  3.6443e-010

参考URL: http://www.doornik.com/support.html#faq_ox

関数で複数の値を返す

関数で複数の値を返すことができます.配列として返すので,波括弧{}で記述し,コンマでそれぞれ区切ります.

#include <oxstd.h>
function(const dX, const dY)
{
  return {dX+dY, dX*dY};
}

main()
{
  decl da, db;
  [da,db]=function(1,2);
  println(da,db);
}

なお,返却値のdbは省略することもできます.

Gaussのプログラムを実行

GaussのプログラムをOx上で実行することができます.拡張子がprgになっていると認識されないので,拡張子をsrcに変更してからOxEditで読み込めば実行できるようになります.

(Uploaded: 2011年06月17日, Last Updated: 2012年08月14日)  ページ上部へ移動

数値微分(多変数関数)

ここでは多変数関数(ƒ:Rn→R)のgradientとHessianを数値的に求めてみます.これは回帰分析では誤差の二乗和Q(β)=u'(β)u(β) (Q:RkR)をパラメータベクトルβに関して微分するのに対応します.

解析的にn階導関数を求めることが難しい場合,コンピュータ上で数値的に導関数を求めることができます.ここでは任意の実数値ベクトルaRnで評価した1階導関数ƒ'(a),2階導関数ƒ''(a)のそれぞれの値を求める方法を紹介します.

簡単化のため,一変数関数で例を示します.もちろん多変数にも拡張可能で,そちらは非線形最小二乗法で扱っています.例として,JavaScriptでニュートン法で用いた関数と同じものを使用します.以下に再掲します.

f(x)

その1階導関数は以下のようになります.

First Derivative of f(x)

2階導関数は以下のようになります.

Second Derivative of f(x)

上記の各導関数をx=3で評価すると,それぞれƒ'(3)=-20,ƒ''(3)=10となります.

Oxにはmaximizeパッケージに1階導関数と2階導関数を求める関数が付属されています.それぞれ,Num1DerivativeNum2Derivativeです.これを用いてx=3における1階導関数,2階導関数の値を数値的に求めてみます.

以下がプログラムです.

#include <oxstd.h>
#import <maximize>

//Num1Derivative, Num2Derivativeで使用する関数を定義
fMyfun( const vP, const adFunc, const avScore, const amHessian ){
  adFunc[0] = vP[0] 4 - 6*vP[0]3 + 5*vP[0]2 + 4*vP[0] + 1;
  return 1;
}

main()
{
  decl dx0 = 3;
  decl df1, df2;
  Num1Derivative( fMyfun, dx0, &df1 );
  Num2Derivative( fMyfun, dx0, &df2 );
  println();
}

プログラムの実行結果は以下のリンクのテキストファイルから閲覧できます.

解析的に1階導関数,2階導関数を求めてから,x=3で評価した値と,数値的に求めた微係数が等しくなっていることがわかります.Num1DerivativeとNum2Derivativeについての詳細はOxのヘルプを確認してください.

【作成したOxのプログラム】: nd_s.ox
【プログラムの実行結果】: output_nd_s.txt

(Uploaded: 2011年02月22日, Last Updated: 2011年02月22日)  ページ上部へ移動

数値微分(多変数ベクトル値関数)

多変数ベクトル値関数(ƒ:RnRm)に対してヤコビ行列を求めたいとします.これは回帰モデルy=x(β)+ux(β)をパラメータベクトルβに関して微分するのに対応します.このような場合,maximizeパッケージのNumJacobianを用います.

ここでは以下の非線形モデルを考えます.

Nonlinear Model

ここで説明変数とパラメータベクトルを以下のスカラー関数として表すことができます.

Scalar Function (Nonlinear Regretion Model).

したがって,多変数ベクトル値関数x(β)のβに関するヤコビ行列の各要素は以下のように計算できます.

Typical Element of Jacobian (beta1),
Typical Element of Jacobian (beta2).

解析的に解いてから評価したヤコビ行列と数値的に求めたヤコビ行列を比較してみます.

ここでx=<1;2;3;4;5>の列ベクトルとし,β=<1;1>で評価したヤコビ行列を数値的に求めます.

以下がプログラムになります.

#include <oxstd.h>
#import <maximize>

//変数をグローバル変数として定義
static decl g_vx;

//NumJacobianで使用する関数を定義
fMyfun( const avR, const vU ){
  avR[0] = vU[0] + vU[1]*exp( - vU[1]*g_vx );
  return 1;
}

main()
{
  //データの作成
  g_vx = <1:5>;

  //NumJacobianのための変数用意
  decl vp0, mjacob;
  decl vp = ones( 2, 1 );

  //NumJacobian
  if( NumJacobian( fMyfun, vp0, &mhess ) ){
    println("Numerical Jacobian at <1; 1> = ", mjacob);
  }
}

プログラムの実行結果は以下のリンクのテキストファイルから閲覧できます.

解析的に求める方は以下のOxのプログラムに付け加えています.比較すると解析的に求めた値と等しくなっていることがわかります.もし解析的に求めることが困難な場合はNumerical Derivativeに頼ることも一つの方法です.

【作成したOxのプログラム】: nd_jacob.ox
【プログラムの実行結果】: output_nd_jacob.txt

(Uploaded: 2011年02月22日, Last Updated: 2011年02月22日)  ページ上部へ移動

非線形最小二乗法 (Newton's Method)

非線形最小二乗法の考え方はOLSと同じで,誤差の二乗和を最小にするパラメータを求めます.しかし,今回は文字通り「非線形」であり,解析的には解くことができません.そこで,Newton's Methodを利用します.

なお,Newton's MethodについてはJavaScriptのページで簡単な解説を載せています(Javascript - JavaScriptでニュートン法).また,ブラウザ上でプログラムを実行できるようにしています.

最小化問題をNewton's Methodを用いて解く時に重要な点は,JavaScriptでニュートン法で述べたように,1変数関数の場合,その2階導関数が常に正であることです.多変数関数の場合では,ここでは回帰分析のケースなので,誤差の2乗和Q(β)=u'(β)u(β) (Q:RkR)のHessianが正値定符号行列になっていることと対応します.

ここでは例として,以下のような非線形回帰モデルを考えてみます.

Nonlinear Regression Model

解析的にgradientとHessianを求めるのが少し困難なケースを選びました.非線形最小二乗法によりパラメータの推定値を求めます.そこで,Newton's Methodを適用します.

JavaScriptでニュートン法の例とは異なり,ここでは多変数関数になっています.したがって,以下の反復式を解くことによって収束する解を求めます.

Iteration for Newton's Method

ここで,H-1(β)はQ(β)のHessianの逆行列,g(β)はQ(β)のgradientを表します.

しかし,収束解が最小値を達成するという保証はなく,初期値をいろいろ試して,その中で最も小さな残差2乗和を達成するパラメータベクトルの推定値を最終的に候補とみなします.

OxにはMaximization用の関数が用意されていますが,経験上,正しい結果が得られたことがあまり多くありません(自分の使い方に問題があるかもしれません).だからと言って,毎回Newton's Methodを使うために,1階導関数,2階導関数を解析的に計算するのは大変です.そこで数値微分(多変数関数)で解説した,Num1Derivative, Num2Derivativeを用い,数値的にgradientとHessianを求めます.

以下がプログラムになります.

#include <oxstd.h>
#include <oxfloat>
#import <maximize>

class MyNewton{
  decl m_vY, m_vX1, m_vX2;
  decl m_cN, m_cPar;
  decl m_asVarNames;
  decl m_vPar;
  decl m_dSSR;
  
  MyNewton();
  LoadVarNames( const asVarNames );
  LoadData( const vY, const mX );
  InitPar( const vPar );
  DoEstimation( vStart );
  Estimate();
  Output();
  fSSR( const vP, const adFunc, const avScore, const amHessian );
};
MyNetwon::MyNewton()
{
  println("Newton's Method \n");
}
MyNetwon::LoadVarNames( const asVarNames )
{
  m_asVarNames = asVarNames;
}
MyNetwon::LoadData( const vY, const mX )
{
  m_cN = rows(vY);
  m_vY = vY;
  m_vX1 = mX[][0];
  m_vX2 = mX[][1];
}
MyNetwon::InitPar( const vPar )
{
  m_vPar = vPar;
  m_cPar = rows(vPar);
}
MyNetwon::DoEstimation( const vStart )
{
  decl vgrd, mhess, mIH;
  decl i, stp, eps, maxIt;
  i = 0;
  eps = 1e-12;
  maxIt = 1e+5;
  
  do{
    if( ++i > maxIt ){
      println("Convergence not achieved.");
      vStart = M_NAN * zeros(m_cPar,1);
      break;
    }
    Num1Derivative( fSSR, vStart, &vgrd );
    Num2Derivative( fSSR, vStart, &mhess );
    mIH = invert( mhess );
    
    vStart = vStart - mIH*vgrd;
    fSSR( vStart, &m_dSSR, 0, 0 );
    
    stp = double( fabs( vgrd'mIH*vgrd ) );
    
    if( isnan(m_dSSR) ){
      vStart = M_NAN * zeros(m_cPar,1);
  }
  println("Iteration ", i, "SSR = ", m_dSSR, "stp = ", stp)
  return vStart;
}
MyNetwon::Estimate()
{
  [m_vPar] = DoEstimation(m_vPar);
  Output();
}
MyNetwon::Output()
{
  println("+++++ Results +++++");
  //以下略
}
MyNewton::fSSR( const vP, const adFunc, const avScore, const amHessian )
{
  decl vu;
  vu = m_vY - vP[0] - (1/vP[1])*log( vP[2]*(m_vX1.^vP[1]) + (1-vP[2])*(m_vX2.^vP[1]) );
  adFunc[0] = double( vu'vu );
  return 1;
}

main()
{
  //データの作成
  ranseed(100);
  decl in, vpt;
  in = 1000;
  vpt = <5; 0.8; 0.4>;
  g_vx1 = 10 + rann( in, 1 );
  g_vx2 = 5*ranu( in, 1 );
  g_vy = vpt[0] + (1/vpt[1])*log( vpt[2]*(g_vx1.^vpt[1]) + (1-vpt[2])*(g_vx2.^vpt[1]) ) + rann( in, 1 );

  decl asVarNames = {"beta1", "beta2", "beta3"};
  decl vp0 = <1; 0.1; 0.1>;

  decl mynewton = new MyNewton();
  mynewton.LoadVarNames( asVarNames );
  mynewton.LoadData( vy, mx );
  mynewton.InitPar( vp0 );
  mynewton.Estimate();
  delete mynewton;
}

プログラムの実行結果は以下のリンクのテキストファイルから閲覧できます.

真の値はvp=<5; 0.8; 0.2>に設定しています.Hessianが常に正値定符号行列でないため,Newton's Methodでは初期値に大きく依存してしまい,最小化を達成する推定値を安定的に求めることができません.事前に真の値に関する情報を全く知らない場合,非常に多くの初期値を試すことになってしまいます.この問題を改善するために次のステップとして非線形最小二乗法 (Levenberg-Marquardt Method)を紹介します.

【作成したOxのプログラム】: nls_newton.ox 
【プログラムの実行結果】: output_nls_newton.txt

(Uploaded: 2011年02月22日, Last Updated: 2011年06月22日)  ページ上部へ移動

非線形最小二乗法 (Gauss-Newton Method)

Hessianが正値定符号を満たさないという問題に対処するために,Gauss-Newton Methodというアルゴリズムがあります.

簡単な説明ですが,Gauss-Newton MethodではHessianをヤコビ行列Jのクロス積J'Jで近似します.理由は下の式ですが,ここではほとんど説明を省略しています.詳細に関して,Davidson and MacKinnon (2004, Chap. 6)の説明が簡潔でわかりやすいと思います.簡単には,βが一致推定量に近づくにつれて,括弧内第一項のyt-xt(β)がゼロに収束していくので,これをゼロとして扱うということです.よってヤコビ行列のクロス積J'Jで近似するということになっています.つまり,括弧内第二項で近似をしています.

Typical Element of Hessian

この場合,正値定符号行列であることは守られますが,極値の周りでの近似になるため,そこから大きく外れた場合は精度が非常に悪いことがわかります.複雑な関数の場合には,なかなか収束しません.

Gauss-Newton Methodの反復計算式は以下になります.

Gauss-Newton Method

上記において,u(β)=y-x(β)は残差ベクトルを表します.ここでは,Newton's Methodと違って符号がプラスに変わっていることに注意してください.

Gauss-Newton Methodを使う別のメリットとして,計算時間の短縮もあげられます.Newton's MethodではHessianを求める必要がありました.推定するパラメータの数が多い場合,数値的にHessianを求めるには非常に時間がかかってしまいます.一方で,Gauss-Newton Methodを用いるためには,残差ベクトルu(β)=y-x(β)とx(β)のβに関するヤコビ行列のみで計算できます.

Gauss-Newton Methodでは収束計算が安定せずうまくいきません.それに対処する方法は非線形最小二乗法 (Levenberg-Marquardt Method)で扱います.

Gauss-Newton Methodのプログラムとその実行結果は以下のリンクのテキストファイルから閲覧できます.

Gauss-Newton Methodの反復計算のJ'(β)u(β)は,gradientにマイナス2分の1倍すれば求めることができます.Gauss-Newton Methodに関してはDavidson and MacKinnon (2004, Chap. 6)に説明があります.

【作成したOxのプログラム】: nls_gn.ox
【プログラムの実行結果】: output_nls_gn.txt

参考文献

  • Davidson, Russell and James G. MacKinnon (2004) Econometric Theory and Methods, New York: Oxford University Press.

(Uploaded: 2011年04月15日, Last Updated: 2011年06月22日)  ページ上部へ移動

非線形最小二乗法 (Levenberg-Marquardt Method)

Gauss-Newton Methodを用いた場合,収束には大きな問題があります.これはHessianの近似方法に問題があります.Hessianに修正を加えたModified Gauss-Newton MethodのことをLevenberg-Marquardt Methodと呼びます.

Levenberg-Marquardt Methodの詳細はChong and Żak (2001), Nocedal and Wright (2006)を参照してください.

Gauss-Newton Methodにシンプルな修正を加えるだけで,かなり安定したアルゴリズムになります.Levenberg-Marquardt Methodの反復計算式は以下になります.

Levenberg-Marquardt Method

ここで,Iは単位行列,cは調整パラメータを表しています.Newton's Methodと違って符号がプラスに変わっていることに注意してください.

ヤコビ行列の積J'Jに単位行列Iを足すだけです.調整パラメータcは収束解に近づくにつれてゼロに近づけていきます.

Levenberg-Marquardt Methodのプログラムとその実行結果は以下のリンクのテキストファイルから閲覧できます.

Hessianの近似であるmjacob'mjacobに対して単位行列unit()を足して修正を行っています.収束解に近づくほどdcの調整パラメータの値をゼロに近づけていきます.

上記ではパラメータベクトルの真の値をvp=<5; 0.8; 0.4>として設定しています.非線形最小二乗法 (Newton's Method)と同じ初期値を試していますが,Levenberg-Marqurdt Methodの方では正しく推定値が得られています.Hessianが正値定符号行列でない場合,Newton's Methodより安定したアルゴリズムになっています.

なお,推定値の共分散行列まで計算していません.Gauss-Newton Regressionという手法を用いて計算できます.また数値微分(多変数ベクトル値関数)で用いたNumJacobianを使うとヤコビ行列を簡単に計算できます.詳細はDavidson and MacKinnon (2004, Chap. 6)を参照してください.

【作成したOxのプログラム】: nls_lm.ox
【プログラムの実行結果】: output_nls_lm.txt

参考文献

  • Chong, Edwin K. and Stanislaw H. Żak (2001) An Introduction to Optimization: New York: John Wiley & Sons, 2nd Ed.
  • Davidson, Russell and James G. MacKinnon (2004) Econometric Theory and Methods, New York: Oxford University Press.
  • Nocedal, Jorge and Stephen J. Wright (2006) Numerical Optimization: New York: Springer, 2nd Ed.

(Uploaded: 2011年04月15日, Last Updated: 2011年06月22日)  ページ上部へ移動

非線形最小二乗法 (Levenberg-Marquardt Modification for Newton's Method)

Levenberg-Marquardt Methodによる修正をそのままNewton's Methodに拡張できます.Hessianを近似しませんが,正値定符号になるように対角要素が全て正の値を持つ対角行列を足します.

Levenberg-Marquardt Modificationの詳細はChong and Żak (2001)を参照してください.

ここでは簡単な説明をします.Newton's MethodでHessianが正値定符号行列をでない場合に探索方向が下方に向かっていきません.最小化とは逆の方向に飛んでしまいます.そこで,以下のようにLevenberg-Marquardt Method流の修正をNetwon's Methodに加えます.

Levenberg-Marquardt Modification for Newton's Method

ここで,Dは対角要素が全て正の値を持つ対角行列,cは調整パラメータを表しています.

Hessianが正値定符号行列になるように対角行列Dを足して修正を加えます.収束する解に近づくにつれて調整パラメータcをゼロに近づけていけば,もとのNetwon's Methodと同じになります.しかし,依然として初期値問題は残ります.いろいろな初期値を試して,最も小さな残差2乗和を達成するパラメータベクトルの推定値を得る必要があります.

Newton's Methodで作成したプログラムにLevenberg-Marquardt Method流の修正を加えています.対角行列Dには単位行列で構いませんが,今回は一様乱数を発生させてHessianの対角要素に足しています.

以下のリンクより,プログラムとその実行結果を閲覧できます.

【作成したOxのプログラム】: nls_lmn.ox
【プログラムの実行結果】: output_nls_lmn.txt

参考文献

  • Chong, Edwin K. and Stanislaw H. Żak (2001) An Introduction to Optimization: New York: John Wiley & Sons, 2nd Ed.

(Uploaded: 2011年02月22日, Last Updated: 2011年06月22日)  ページ上部へ移動

離散選択モデル

離散選択モデルの推定方法.

プロビット

ロジット

順序ロジット

(Uploaded: 2011年06月09日, Last Updated: 2011年06月09日)  ページ上部へ移動

プログラミングのTips

クラス

継承と仮想関数

(Uploaded: 2011年06月22日, Last Updated: 2011年06月22日)  ページ上部へ移動

数値解析(最適化)

  • Chong, Edwin K. and Stanislaw H. Żak (2001) An Introduction to Optimization: New York: John Wiley & Sons, 2nd Ed.
  • Nocedal, Jorge and Stephen J. Wright (2006) Numerical Optimization: New York: Springer, 2nd Ed.
  • 金谷健一 (2005)『これなら分かる最適化数学:基礎原理から計算手法まで』,共立出版,東京.

線形モデルと非線形モデルの間には大きな障壁があり,非線形モデルになると急激に難しくなります.厳密には理論を理解できていませんので,そのような状態で実用しようとすると計算手法上で行き詰ったり,何らかの問題に遭遇します.そんな時,上記の本の説明がとても参考になっています.

Chong and Żak (2001),Nocedal and Wright (2006)は詳しく書かれていると思います.金谷 (2005)はコンパクトにまとまっています.

プログラミング言語

  • Doornik, Jurgen A. and Marius Ooms (2006) Introduction to Ox: http://www.doornik.com/
  • 柴田望洋 (2004) 『明解C言語 入門編』,ソフトバンククリエイティブ,東京,新版.
  • 柴田望洋 (2004) 『明解C言語 中級編』,ソフトバンククリエイティブ,東京,新版.
  • 柴田望洋 (2007) 『明解Java 入門編』,ソフトバンククリエイティブ,東京,新版.
  • 柴田望洋 (2009) 『明解C++ 入門編』,ソフトバンククリエイティブ,東京,新版.
  • 矢沢久雄 (2002) 『C++ クラスと継承完全制覇』,技術評論社 ,東京.

C/C++やJavaとほぼ同じ構造なので上記の本で学んだらOxへの理解も深まると思います.

(Uploaded: 2011年02月22日, Last Updated: 2016年01月20日)  ページ上部へ移動