2006-12-12

累積確率密度分布関数の近似

12月11日の記事で、線形近似式が出てきて、その式が、多項展開とその積分の便宜のために、ちょっと複雑になっていた。

また、累積密度関数の近似式を得る作業の途中であり、(0,0),(1,1)を通るような線形近似式を得るためのソースが必要になっている。

準備１
- n次多項式による最小自乗法では、(n+1)個の係数 $a_i$ , $i=0,1,...,n$ によって表された関数 $F(x;a_0,a_1,...,a_n)=￥sum_{i=0}^{n}a_ix^i$ がある。
- 今、,,なるデータがあったとき、最小化したい値は
  - $J=￥sum_{u=1}^{U}(y_u-F(x_u))^2$ である。
- このを最小化するようなは、それぞれの偏微分が0となるから
  - なるn+1個の等式が得られる
    - $￥frac{￥partial J}{￥partial a_i}=-2￥sum_{u=1}^{U}(y_u-F(x_u))￥frac{￥partial F(x_u)}{￥partial a_i}=0$ , $￥frac{￥partial F(x_u)}{￥partial a_i}=x_u^i$
- このn+1個のに関する連立1次方程式を解くことで、が得られる。一般に、k変数の連立1次方程式はkxkの行列Mと長さkのベクトルA,Bとでと表記できる。ただし、Aは変数に相当する。上記の偏微分に基づく連立方程式の場合には
  - $A=(a_0,a_1,...,a_n)$ であり、
  - $M=(m_{i,j})$ , $m_{i,j}=￥sum_{u=1}^{U}x_u^jx_u^i$ $i,j=0,1,...,n$ 、
  - $B=(b_i)$ , $b_i=￥sum_{u=1}^{U}x_u^iy_u$ , $i=0,1,...,n$
- の逆行列を求めることで、としてを求められることがわかるから、計算ツールとして必要なのは、任意の次数nに対して、データから行列Mを作るサブルーチンと、逆行列を作るサブルーチンと、行列の積を求めるサブルーチンの３つであることがわかる。
  - 行列の積(ベクターも行列扱いしている)



	public static double[][] Prodmat(double[][] a, double[][] b){

		double[][] ret = new double[a.length][b[0].length];

		for(int i=0;i<ret.length;i++){

			for(int j=0;j<ret[i].length;j++){

				ret[i][j]=0;

				for(int k=0;k<a[i].length;k++){

					ret[i][j]+=a[i][k]*b[k][j];

				}

			}

		}

		return ret;

	}

逆行列のソースはこちら
次数とデータセットを引数にし、行列Mを作り、その逆行列を作成した上で $A=M^{-1}B$ を行うことで、n次関数のn+1個の係数を返すソース



	public static double[] LstSq(int n,double[] x, double[] y){

//n次式、x,yはデータセット

		double[] ret = new double[n+1];

		double[][] a = new double[n+1][n+1];

		double[][] Y = new double[n+1][1];

		

		double[][] xK = new double[x.length][n+1];

		for(int i=0;i<xK.length;i++){

			xK[i][0] = 1;

			for(int j=1;j<xK[i].length;j++){

				xK[i][j] = xK[i][j-1]*x[i];

				//System.out.println("i j " + i + " " + j +" " + xK[i][j]);

			}

		}

		

		

		for(int i=0;i<a.length;i++){

			for(int j=0;j<a[i].length;j++){

				for(int k=0;k<x.length;k++){

					a[i][j] +=xK[k][i]*xK[k][j];

				}

			}

		}

		for(int i=0;i<Y.length;i++){

			for(int j=0;j<x.length;j++){

				Y[i][0] += y[j]*xK[j][i];

			}

		}

		

		double[][] X = Renritsu(a,Y);

		//System.out.println("X.leng " + X.length);

		for(int i=0;i<X.length;i++){

			for(int j=0;j<X[i].length;j++){

				//System.out.println(X[i][j]);

			}

		}

		for(int i=0;i<ret.length;i++){

			ret[i] = X[i][0];

		}

		//ret = X[0];

		return ret;

	}

累積確率密度分布であるための条件を取り込む
- 累積確率密度分布関数は、以下の条件を持つ
  - (1) 0から1の間で定義されること
  - (2) (0,0)を通る
  - (3) (1,1)を通ること
  - (4) 非減少関数であること
- この４条件のうち、
  - (1)は定義範囲での近似がよければ、その外側の挙動には関知しなくてよい、ということで近似作業にあたって無視することとする。
  - (4)は、観測データの与え方を工夫することによって、特に配慮しなくとも、実現されることを期待する。具体的には、十分に密なデータポイントを与え、しかも、そのデータポイントが非減少であるようにすることで、その期待は十分に現実的であると予想する。
  - (2)は、 $F(x)=￥sum_{i=1}^{n}a_i x^i$ のように、0次の項をなくす(n次多項式にて係数n個、偏微分式n個、nxn行列の作成：n次連立方程式)ことで対応できる
  - (3)は、という条件を持ち込む。この条件よりと、n次の項の係数が、の関数となっていることがわかる。このことから、n-1次の連立方程式となり、
    - と変化し、そのことから、(n-1)x(n-1)行列Mと、長さnのベクトルBとは、それぞれ
      - $M=(m_{i,j})$ , $m_{i,j}=￥sum_{u=1}^{U}(x_u^j-x_u^n)(x_u^i-x_u^n)$ , $i,j=1,2,...n-1$
      - $B=(b_i)$ , $b_i=￥sum_{u=1}^{U}(x_u^i-x_u^n)y_u$ , $i=1,2,...n-1$



public static double[] LstSq0011_2(int n,double[] x, double[] y){

		n--;//変数を1つ減らす

		double[] ret = new double[n+2];

		double[][] a = new double[n][n];

		double[][] Y = new double[n][1];

		

		double[][] xK = new double[x.length][n+1];

		for(int i=0;i<xK.length;i++){

			//xK[i][0] = 1;

			xK[i][0] = x[i];

			for(int j=1;j<xK[i].length;j++){

				xK[i][j] = xK[i][j-1]*x[i];

				//System.out.println("i j " + i + " " + j +" " + xK[i][j]);

			}

		}

		

		

		for(int i=0;i<a.length;i++){

			for(int j=0;j<a[i].length;j++){

				for(int k=0;k<x.length;k++){

					a[i][j] +=(xK[k][i]-xK[k][n])*(xK[k][j]-xK[k][n]);

				}

			}

		}

		//for(int i=0;i<a[a.length-1].length;i++){

			//a[a.length-1][i]=1;

		//}

		for(int i=0;i<Y.length;i++){

			for(int j=0;j<x.length;j++){

				Y[i][0] += (y[j]-xK[j][n])*(xK[j][i]-xK[j][n]);

			}

		}

		//Y[Y.length-1][0]=1;

		double[][] X = Renritsu(a,Y);

		//System.out.println("X.leng " + X.length);

		for(int i=0;i<X.length;i++){

			for(int j=0;j<X[i].length;j++){

				//System.out.println(X[i][j]);

			}

		}

		ret[0]=0;

		ret[ret.length-1]=1;

		for(int i=1;i<ret.length-1;i++){

			ret[i] = X[i-1][0];

			ret[ret.length-1]-=ret[i];

		}

		

		

		//ret = X[0];

		return ret;

	}

上述のようにして $F(x)=￥sum_{i=1}^{n}a_ix^i$ の形で得られた式を、 $F(x)=1-￥sum_{k=1}^{K} c_k (1-x)^k$ という式に変形するにあたっては、-- $a_k=(-1)^{k+1}￥sum_{i=k}^{K}(c_i ￥times _iC_k)$ なる関係にあることを昨日の記事で示してあるので、この関係を $A=(a_i),C=(c_i)$ として $A=MC$ という連立方程式とすることで、上記、最小自乗法のときと同様に解くことで得られる。



public static double[] StructureAreaTrans(double[] d){

		int len = d.length;

		double[][] a = new double[len][len];

		double[][] dm = new double[len][1];

		for(int i =0;i<d.length;i++){

			dm[i][0] = d[i];

		}

		//double pm = -1;

		a[0][0] = 1;

		for(int i=1;i<a[0].length;i++){

			a[0][i] = a[0][i-1]*(double)(i+1)/(double)(i+1-1);

			//System.out.println("i j a[i][j] " + "0 " + i + " " + a[0][i]);

		}

		for(int i=1;i<a.length;i++){

			//pm *= -1;

			for(int j=0;j<a[i].length;j++){

				if(j<i){

					a[i][j]=0;

				}else{

					a[i][j] = (-1)*a[i-1][j-1]*(double)(j+1)/(double)(i+1);

				}

				//System.out.println("i j a[i][j] "+ i +" " + j + " " + a[i][j]);

			}

		}

		double[][] ag = GyakuGyouretsu.inv(a);

		double[][] am = Prodmat(ag,dm);

		double[] ret = new double[len];

		for(int i=0;i<ret.length;i++){

			ret[i] = am[i][0];

		}

		return ret;

	}