/*  -*- MaTX -*-
 *
 *  【名前】
 *      hist() - ヒストグラム
 *
 *  【形式】
 *      {N,X} = hist(Y)
 *         Array N,X;
 *         Array Y;
 *
 *      {N,X} = hist(Y,k)
 *         Array N,X;
 *         Array Y;
 *         Integer k;
 *
 *  【機能説明】
 *      {N,X} = hist(Y)は，Y の全てのデータの最小値と最大値の間を 10
 *      等分し，そのヒストグラムを表示するために必要な頻度 N とその位置
 *      を示すベクトル X を返す。得られた X と N を使って，barp(X,N) で
 *      ヒストグラムをプロットできる。hist(Y,k) は最小値と最大値の間を
 *      k 等分する。
 * 
 *  【例題】
 *      正規分布するデータから鐘状のヒストグラムを作る。
 *
 *      Y = randn(1,1000);
 *      {N,X} = hist(Y,40);
 *      barp(X,N)
 *
 *  【関連項目】
 *      hist_row(), hist_col(), bar(), barp()
 *
 *   ------------------------------------------------------------------------
 *
 *  【名前】
 *      hist_col() - 列毎のヒストグラム
 *
 *  【形式】
 *      {N,X} = hist_col(Y)
 *         Array N,X;
 *         Array Y;
 *
 *      {N,X} = hist_col(Y,k)
 *         Array N,X;
 *         Array Y;
 *         Integer k;
 *
 *  【機能説明】
 *      {N,X} = hist_col(Y)は，ベクトル Y データの最小値と最大値の間を 10
 *      等分し，そのヒストグラムを表示するために必要な頻度 N とその位置
 *      を示すベクトル X を返す。得られた X と N を使って，barp(X,N) で
 *      ヒストグラムをプロットできる。hist_col(Y,k) は最小値と最大値の間を
 *      k 等分する。
 * 
 *      Y が行列のとき，列毎にヒストグラムを表示するために必要な頻度 N と
 *      とその位置を示す X を返す。行列 X と N の対応する列を使って，
 *      barp(X(:,i),N(:,i)) でヒストグラムをプロットできる。
 *
 *  【例題】
 *      正規分布するデータから鐘状のヒストグラムを作る。
 *
 *      Y = randn(1000,2);
 *      {N,X} = hist_col(Y,40);
 *      barp(X(:,1),N(:,1));
 *
 *  【関連項目】
 *      hist(), hist_row(), bar(), barp()
 *
 *   ------------------------------------------------------------------------
 *
 *  【名前】
 *      hist_row() - 行毎のヒストグラム
 *
 *  【形式】
 *      {N,X} = hist_row(Y)
 *         Array N,X;
 *         Array Y;
 *
 *      {N,X} = hist_row(Y,k)
 *         Array N,X;
 *         Array Y;
 *         Integer k;
 *
 *  【機能説明】
 *      {N,X} = hist_row(Y)は，ベクトル Y データの最小値と最大値の間を 10
 *      等分し，そのヒストグラムを表示するために必要な頻度 N とその位置
 *      を示すベクトル X を返す。得られた X と N を使って，barp(X,N) で
 *      ヒストグラムをプロットできる。hist_row(Y,k) は最小値と最大値の間を
 *      k 等分する。
 * 
 *      Y が行列のとき，行毎にヒストグラムを表示するために必要な頻度 N と
 *      とその位置を示す X を返す。行列 X と N の対応する行を使って，
 *      barp(X(i,:),N(i,:)) でヒストグラムをプロットできる。
 *
 *  【例題】
 *      正規分布するデータから鐘状のヒストグラムを作る。
 *
 *      Y = randn(2,1000);
 *      {N,X} = hist_row(Y,40);
 *      barp(X(1,:),N(1,:));
 *
 *  【関連項目】
 *      hist(), hist_col(), bar(), barp()
 */

Func List hist(Y, n, ...)
	Array Y;
	Integer n;
{
	error(nargchk(1, 2, nargs, "hist"));
	
	if (nargs == 1) {
		return hist_row(reshape(Y, 1, Rows(Y)*Cols(Y)));
	} else {
		return hist_row(reshape(Y, 1, Rows(Y)*Cols(Y)), n);
	}
}

Func List hist_col(Y_, n, ...)
	Array Y_;
	Integer n;
{
	Integer i,j;
	Array y,nn,x,xx,N,X,Y;
	Real y_max, y_min, bw;
	
	error(nargchk(1, 2, nargs, "hist_col"));

	if (nargs == 1) {
		n = 10;
	}
	
	if (Rows(Y_) == 1 && Cols(Y_) != 1) {
		Y = makecolv(Y_);
	} else {
		Y = Y_;
	}

	X = N = Z(n,Cols(Y));
	for (i = 1; i <= Cols(Y); i++) {
		y = Y(:,i);

		y_min = min(y);
		y_max = max(y);
		bw = (y_max - y_min)/n;
		x = trans(y_min .+ bw*[0:n]);
		x(length(x)) = y_max;

		nn = Z(n+1,1);
		for (j = 2; j <= n + 1; j++) {
			nn(j,1) = length(find(y <= x(j)));
		}

		N(:,i) = nn(2:n+1,1) - nn(1:n,1);
		X(:,i) = x(1:length(x)-1) .+ bw/2;
	}

	if (Rows(Y_) == 1 && Cols(Y_) != 1) {
		N = trans(N);
		X = trans(X);
	}

	return {N, X};
}

Func List hist_row(Y, n, ...)
	Array Y;
	Integer n;
{
	Array N, X;

	error(nargchk(1, 2, nargs, "hist_row"));
	
	if (nargs == 1) {
		{N, X} = hist_col(trans(Matrix(Y)));
	} else {
		{N, X} = hist_col(trans(Matrix(Y)), n);
	}
	N = trans(Matrix(N));
	X = trans(Matrix(X));
	return {N, X};
}