Classiclll's Blog

an old boy

ロジスティックモデル同定　ベータテスト８の解説と終了宣言

クラス遊び #その他趣味

VfuncとEqSysとMatとVecの例題： ロジスティックモデル同定



妄想中のVec, Mat、VfuncとEqSysのベータテスト最終ケースにロジスティックモデルの同定を選んだ。

ロジスティックモデルは、基本的な成長モデルであり、例えば人口予測の数学モデル等に応用されている。
　　取り上げたモデル式は、y = a / (b + c^x)

【予測】の一般的仮定
　世の中の（数量で捉えられる）気になる事象は、何らかの数式で規定される。しかし、我々がその事象を観測する時、様々な原因から観測誤差が混ざった形でしか計測できない。多くの観測を集めれば、観測誤差をある程度取り除く事が出来るはず。
　　y：観測値　＝　f(x)：規定する数式　＋　ε：観測誤差
　　　＊x：観測値を得たときの前提（説明変量）
　　　＊ε：観測誤差の期待値はゼロ

【ロジスティックモデル】への期待
　個体が少ない社会では生殖の機会が少なく、個体数の増加率は少ない。その後、若い個体の増加に従い生殖の機会が増加し、個体数は急激に増加する。しかし、個体数の急激な増加は居住面積を含めた様々な個体存在に必須なリソースが消費され、生殖は抑制され、個体数の増加が鈍くなり、最後には出生数と死亡数が拮抗し、個体数の増加は止まる。
　このような考察に良く合う曲線がロジスティックなので、もしかして人口増加の推移を良く説明してくれるのではないか？

【観測に基づくパラメータ推定】
　多数の観測値(obs)と観測したときの説明変量(samples)が得られているとする。このとき全ての観測誤差の累計を最小とするパラメータ(a,b,c)の組を得たい。
なので、例えば
　　ssr(a,b,c) = Σ(obs[i]-f(samples[i])^2 ：　残差平方和（SSR)
を最小にする(a,b,c)が所望のものとなる。

。α、β、γが所望のものとすると、下記を満たす。（必要条件）
　　{ssr(α,β,γ)}/d(a,b,c)=0　⇔　Σ2*(f(samples[j])-obs[j])*{f(samples[j])}/d(a,b,c) = 0
クラスVfunは上記の定式化は実装済みなので
　　{ f(x) ]/da = 0　　⇔　 1 / (β+γ^x) = 0
　　{ f(x) ]/db = 0　　⇔　 -α / (β+γ^x)^2 = 0
　　{ f(x) ]/dc = 0　　⇔　 -α*x*γ^(x-1) / (β+γ^x)^2 = 0
を継承クラスModifiedLogisticに実装し、上記非線形連立方程式をVfuncに解かせればOK
ここまでが、ケース８のためのテストドライバと実装クラスの説明

テスト結果を見ると、ちょっと無理っぽい
　・Simplexは収束しない
　・Newtonはウソっぽい（実際、上記式では途中でダメになる）

ま、テスト５で例のあれが解ける事が確認できたので、とりあえずここまで。

==== ベータテストドライバ （８）====
//
//  Test for Vfunc, EqSys, Mat, Vec, and their solver methods.
//  Created by Classiclll on 06/03/19.
//

import lll.Loc.*;
public void setup() {
//	 8. estimate multi-variated non-linear regression coeficient, 
//		f8(x) = a / (b + c^x) - modified Logistic curve 
  Vec obs = new Vec(100);
  Mat samples = new Mat(100,2);
  ModifiedLogistic f8 = new ModifiedLogistic(); // f8(x) = a / (b + c^x)
  Vec trueParam = new Vec(new double[]{3,5,1});
  for (int i=0;i<100;i++) { // generate the random samples
	  double[] smpl 
	  	= { (Math.random()-0.5)*5, (Math.random()-0.5)*5};
	  samples.setRowVec(new Vec(smpl),i);
	  obs.arrayRef()[i]
	     = f8.valueAt(samples.rowVec(i), trueParam) + (Math.random()-0.5)*1;
  }
  Vec p0 = new Vec(new double[]{5,5,5});
 
  println("\n[ multi-variated non-linear regression Test ]");
  println(">> f8(1), value of f1 at x=1,y=1,a=3,b=5,c=1 is "
		  + f8.valueAt(new Vec(new double[]{1}), trueParam));
  println(">>true param is" + trueParam);
  println(">>estimated param by Newton "+f8.bestPrmByNewton(p0, obs, samples));
  println(">>estimated param by Simplex "+f8.bestPrmBySimplex(p0, obs, samples, 5000));

テスト結果

[ multi-variated non-linear regression Test ]
>> f8(1), value of f1 at x=1,,a=3,b=5,c=1 is 0.5
>>true param isVec(3.0 ,5.0 ,1.0)
>>estimated param by Newton Vec(3.0337683405566245 ,5.0 ,1.1491032007315607)
>>estimated param by Simplex Vec(NaN)

Vfunc、EqSysのテスト用実装クラス

import lll.Loc.*;

//8. estimate multi-variated non-linear regression coeficient, 
//f8(x) = a / (b + c^x) - modified Logistic curve 

public class ModifiedLogistic extends Vfunc {
  public ModifiedLogistic() {
	super(1,3);// dim of domain is 1, dim of parameters is 3
  }
  public double valueAt(Vec x, Vec p) {	// f(x) = a / (b + c^x) 
	return p.elem(0) / (p.elem(1) + Math.pow(p.elem(2),x.elem(0)) );
  }
  public Vec gradAt(Vec x, Vec p) {	// f'(x) = (-a*ln(c)/(b+c^x)^2)
	return this.diffAt(x, p.mul(.01), p);
  }
  public Vec gradParamAt(Vec x, Vec p) {	
	// fp'(p) = (1,   -1/b^2,   -x*c^(x-1)/c^2x )
	double xx=x.elem(0), b=p.elem(1), c=p.elem(2);
	  return new Vec(
	    new double[]{1,  -1/Math.pow(b, 2),
    			-xx*Math.pow(c, xx-1) / Math.pow(c, 2*xx) } );
  }
/*
  public Vec gradParamAt(Vec x, Vec p) {	
	// f(p) = f(a,b,c,d) = a / (b + c^x)
	// fp'(p) = ( 1/(b+c^x),
	//			-a/(b+c^x)^2,
	//			-a*x*c^(x-1)/(b+c^x+d^y)^2  )
	double xx=x.elem(0), a=p.elem(0), b=p.elem(1), c=p.elem(2);
	double denomi = b + Math.pow(c, xx);
	return new Vec(
		new double[]{1/denomi, -a/denomi/denomi,
		-a*xx/denomi*Math.pow(c, xx-1)/denomi } );
  }
*/
}