静岡理工科大学 菅沼ホーム 目次 索引

最小二乗法

    1. A. C++
    2. B. Java
    3. C. JavaScript
    4. D. PHP
    5. E. Ruby
    6. F. Python
    7. G. C#
    8. H. VB

  プログラムは,最小二乗法(多項式近似)を実行するためのものです.データの入力方法に関しては,プログラムの最後に提示してある入力例に対する説明を参考にして下さい.

  1. C++

    /****************************/
    /* 最小二乗法(多項式近似) */
    /*      coded by Y.Suganuma */
    /****************************/
    #include <stdio.h>
    #include <math.h>
    
    int gauss(double **, int, int, double);
    double *least(int, int, double *, double *);
    
    int main()
    {
    	double *x, *y, *z;
    	int i1, m, n;
    
    	scanf("%d %d", &m, &n);   // 多項式の次数とデータの数
    
    	x = new double [n];
    	y = new double [n];
    
    	for (i1 = 0; i1 < n; i1++)   // データ
    		scanf("%lf %lf", &x[i1], &y[i1]);
    
    	z = least(m, n, x, y);
    
    	if (z != NULL) {
    		printf("結果\n");
    		for (i1 = 0; i1 < m+1; i1++)
    			printf("   %d 次の係数 %f\n", m-i1, z[i1]);
    		delete [] z;
    	}
    	else
    		printf("***error  逆行列を求めることができませんでした\n");
    
    	delete [] x;
    	delete [] y;
    
    	return 0;
    }
    
    /******************************************/
    /* 最初2乗法                             */
    /*      m : 多項式の次数                  */
    /*      n : データの数                    */
    /*      x,y : データ                      */
    /*      return : 多項式の係数(高次から) */
    /*               エラーの場合はNULLを返す */
    /******************************************/
    double *least(int m, int n, double *x, double *y)
    {
    	double **A, **w, *z, x1, x2;
    	int i1, i2, i3, sw;
    
    	m++;
    	z = new double [m];
    	w = new double * [m];
    	for (i1 = 0; i1 < m; i1++)
    		w[i1] = new double [m+1];
    	A = new double * [n];
    
    	for (i1 = 0; i1 < n; i1++) {
    		A[i1]      = new double [m];
    		A[i1][m-2] = x[i1];
    		A[i1][m-1] = 1.0;
    		x1         = A[i1][m-2];
    		x2         = x1;
    		for (i2 = m-3; i2 >= 0; i2--) {
    			x2        *= x1;
    			A[i1][i2]  = x2;
    		}
    	}
    
    	for (i1 = 0; i1 < m; i1++) {
    		for (i2 = 0; i2 < m; i2++) {
    			w[i1][i2] = 0.0;
    			for (i3 = 0; i3 < n; i3++)
    				w[i1][i2] += A[i3][i1] * A[i3][i2];
    		}
    	}
    
    	for (i1 = 0; i1 < m; i1++) {
    		w[i1][m] = 0.0;
    		for (i2 = 0; i2 < n; i2++)
    			w[i1][m] += A[i2][i1] * y[i2];
    	}
    
    	sw = gauss(w, m, 1, 1.0e-10);
    
    	if (sw == 0) {
    		for (i1 = 0; i1 < m; i1++)
    			z[i1] = w[i1][m];
    	}
    	else
    		z = NULL;
    
    	for (i1 = 0; i1 < n; i1++)
    		delete [] A[i1];
    	for (i1 = 0; i1 < m; i1++)
    		delete [] w[i1];
    	delete [] A;
    	delete [] w;
    
    	return z;
    }
    
    /*******************************************************/
    /* 線形連立方程式を解く(逆行列を求める)              */
    /*      w : 方程式の左辺及び右辺                       */
    /*      n : 方程式の数                                 */
    /*      m : 方程式の右辺の列の数                       */
    /*      eps : 正則性を判定する規準                     */
    /*      return : =0 : 正常                             */
    /*               =1 : 逆行列が存在しない               */
    /*******************************************************/
    int gauss(double **w, int n, int m, double eps)
    {
    	double y1, y2;
    	int ind = 0, nm, m1, m2, i1, i2, i3;
    
    	nm = n + m;
    
    	for (i1 = 0; i1 < n && ind == 0; i1++) {
    
    		y1 = .0;
    		m1 = i1 + 1;
    		m2 = 0;
    
    		for (i2 = i1; i2 < n; i2++) {
    			y2 = fabs(w[i2][i1]);
    			if (y1 < y2) {
    				y1 = y2;
    				m2 = i2;
    			}
    		}
    
    		if (y1 < eps)
    			ind = 1;
    
    		else {
    
    			for (i2 = i1; i2 < nm; i2++) {
    				y1        = w[i1][i2];
    				w[i1][i2] = w[m2][i2];
    				w[m2][i2] = y1;
    			}
    
    			y1 = 1.0 / w[i1][i1];
    
    			for (i2 = m1; i2 < nm; i2++)
    				w[i1][i2] *= y1;
    
    			for (i2 = 0; i2 < n; i2++) {
    				if (i2 != i1) {
    					for (i3 = m1; i3 < nm; i3++)
    						w[i2][i3] -= w[i2][i1] * w[i1][i3];
    				}
    			}
    		}
    	}
    
    	return(ind);
    }
    
    /*
    -----------データ例1(コメント部分を除いて下さい)---------
    1 10   // 多項式の次数とデータの数
    0  2.450   // 以下,(x, y)
    1  2.615
    2  3.276
    3  3.294
    4  3.778
    5  4.009
    6  3.920
    7  4.267
    8  4.805
    9  5.656
    
    -------------------------データ例2-------------------------
    2 6
    10  28.2
    15  47.0
    20  44.4
    26  32.8
    32  20.8
    40   0.8
    */
    			

  2. Java

      アプレット版では,任意のデータに対して画面上で実行することができます.
    /****************************/
    /* 最小二乗法(多項式近似) */
    /*      coded by Y.Suganuma */
    /****************************/
    import java.io.*;
    import java.util.StringTokenizer;
    
    public class Test {
    	public static void main(String args[]) throws IOException
    	{
    		double x[], y[], z[];
    		int i1, m, n, sw;
    		StringTokenizer str;
    		BufferedReader in = new BufferedReader(new InputStreamReader(System.in));
    					// 多項式の次数とデータの数
    		str = new StringTokenizer(in.readLine(), " ");
    		m   = Integer.parseInt(str.nextToken());
    		n   = Integer.parseInt(str.nextToken());
    
    		x = new double [n];
    		y = new double [n];
    		z = new double [m+1];
    					// データ
    		for (i1 = 0; i1 < n; i1++) {
    			str = new StringTokenizer(in.readLine(), " ");
    			x[i1] = Double.parseDouble(str.nextToken());
    			y[i1] = Double.parseDouble(str.nextToken());
    		}
    
    		sw = least(m, n, x, y, z);
    
    		if (sw == 0) {
    			System.out.println("結果");
    			for (i1 = 0; i1 < m+1; i1++)
    				System.out.println("   " + (m-i1) + " 次の係数 " + z[i1]);
    		}
    		else
    			System.out.println("***error  逆行列を求めることができませんでした");
    	}
    
    	/*************************************/
    	/* 最初2乗法                        */
    	/*      m : 多項式の次数             */
    	/*      n : データの数               */
    	/*      x,y : データ                 */
    	/*      z : 多項式の係数(高次から) */
    	/*      return : =0 : 正常           */
    	/*               =1 : エラー         */
    	/*************************************/
    	static int least(int m, int n, double x[], double y[], double z[])
    	{
    		double A[][], w[][], x1, x2;
    		int i1, i2, i3, sw = 0;
    
    		m++;
    		w = new double [m][m+1];
    		A = new double [n][m];
    
    		for (i1 = 0; i1 < n; i1++) {
    			A[i1][m-2] = x[i1];
    			A[i1][m-1] = 1.0;
    			x1         = A[i1][m-2];
    			x2         = x1;
    			for (i2 = m-3; i2 >= 0; i2--) {
    				x2        *= x1;
    				A[i1][i2]  = x2;
    			}
    		}
    
    		for (i1 = 0; i1 < m; i1++) {
    			for (i2 = 0; i2 < m; i2++) {
    				w[i1][i2] = 0.0;
    				for (i3 = 0; i3 < n; i3++)
    					w[i1][i2] += A[i3][i1] * A[i3][i2];
    			}
    		}
    
    		for (i1 = 0; i1 < m; i1++) {
    			w[i1][m] = 0.0;
    			for (i2 = 0; i2 < n; i2++)
    				w[i1][m] += A[i2][i1] * y[i2];
    		}
    
    		sw = gauss(w, m, 1, 1.0e-10);
    
    		if (sw == 0) {
    			for (i1 = 0; i1 < m; i1++)
    				z[i1] = w[i1][m];
    		}
    		else
    			sw = 1;
    
    		return sw;
    	}
    
    	/*******************************************************/
    	/* 線形連立方程式を解く(逆行列を求める)              */
    	/*      w : 方程式の左辺及び右辺                       */
    	/*      n : 方程式の数                                 */
    	/*      m : 方程式の右辺の列の数                       */
    	/*      eps : 逆行列の存在を判定する規準               */
    	/*      return : =0 : 正常                             */
    	/*               =1 : 逆行列が存在しない               */
    	/*******************************************************/
    	static int gauss(double w[][], int n, int m, double eps)
    	{
    		double y1, y2;
    		int ind = 0, nm, m1, m2, i1, i2, i3;
    
    		nm = n + m;
    
    		for (i1 = 0; i1 < n && ind == 0; i1++) {
    
    			y1 = .0;
    			m1 = i1 + 1;
    			m2 = 0;
    						// ピボット要素の選択
    			for (i2 = i1; i2 < n; i2++) {
    				y2 = Math.abs(w[i2][i1]);
    				if (y1 < y2) {
    					y1 = y2;
    					m2 = i2;
    				}
    			}
    						// 逆行列が存在しない
    			if (y1 < eps)
    				ind = 1;
    						// 逆行列が存在する
    			else {
    							// 行の入れ替え
    				for (i2 = i1; i2 < nm; i2++) {
    					y1        = w[i1][i2];
    					w[i1][i2] = w[m2][i2];
    					w[m2][i2] = y1;
    				}
    							// 掃き出し操作
    				y1 = 1.0 / w[i1][i1];
    
    				for (i2 = m1; i2 < nm; i2++)
    					w[i1][i2] *= y1;
    
    				for (i2 = 0; i2 < n; i2++) {
    					if (i2 != i1) {
    						for (i3 = m1; i3 < nm; i3++)
    							w[i2][i3] -= w[i2][i1] * w[i1][i3];
    					}
    				}
    			}
    		}
    
    		return(ind);
    	}
    }
    
    /*
    -----------データ例1(コメント部分を除いて下さい)---------
    1 10   // 多項式の次数とデータの数
    0  2.450   // 以下,(x, y)
    1  2.615
    2  3.276
    3  3.294
    4  3.778
    5  4.009
    6  3.920
    7  4.267
    8  4.805
    9  5.656
    
    -------------------------データ例2-------------------------
    2 6
    10  28.2
    15  47.0
    20  44.4
    26  32.8
    32  20.8
    40   0.8
    */
    			
    アプレット版
    /****************************/
    /* 最小二乗法(多項式近似) */
    /*      coded by Y.Suganuma */
    /****************************/
    import java.awt.*;
    import java.awt.event.*;
    import java.applet.*;
    import java.util.StringTokenizer;
    
    public class Least extends Applet implements ActionListener {
    	TextField tx1, tx2;
    	TextArea ta1, ta2;
    	Button bt;
    	Font f = new Font("TimesRoman", Font.BOLD, 20);
    
    	/************/
    	/* 初期設定 */
    	/************/
    	public void init()
    	{
    					// レイアウト,背景色,フォント
    		setLayout(new BorderLayout(5, 5));
    		setBackground(new Color(225, 255, 225));
    		setFont(f);
    					// 上のパネル
    		Panel pn1 = new Panel();
    		add(pn1, BorderLayout.NORTH);
    
    		pn1.add(new Label("多項式の次数:"));
    		tx1 = new TextField("1", 2);
    		pn1.add(tx1);
    
    		pn1.add(new Label("データの数:"));
    		tx2 = new TextField("10", 3);
    		pn1.add(tx2);
    					// 中央のパネル(データ)
    		Panel pn2 = new Panel();
    		add(pn2, BorderLayout.CENTER);
    
    		pn2.add(new Label("データ(x y):"));
    		ta1 = new TextArea("0  2.450\n1  2.615\n2  3.276\n3  3.294\n4  3.778\n5  4.009\n6  3.920\n7  4.267\n8  4.805\n9  5.656\n", 15, 15);
    		pn2.add(ta1);
    
    		pn2.add(new Label(" "));
    		bt = new Button("OK");
    		bt.setBackground(Color.pink);
    		bt.addActionListener(this);
    		pn2.add(bt);
    					// 下のパネル(結果)
    		Panel pn3 = new Panel();
    		add(pn3, BorderLayout.SOUTH);
    
    		pn3.add(new Label("係数(次数の高い順):"));
    		ta2 = new TextArea(7, 25);
    		pn3.add(ta2);
    	}
    
    	/******************************/
    	/* ボタンが押されたときの処理 */
    	/******************************/
    	public void actionPerformed(ActionEvent e)
    	{
    		if (e.getSource() == bt) {
    					// データ
    			int m = Integer.parseInt(tx1.getText());
    			int n = Integer.parseInt(tx2.getText());
    			double x[] = new double [n];
    			double y[] = new double [n];
    			double z[] = new double [m+1];
    
    			String s = ta1.getText();
    			StringTokenizer str = new StringTokenizer(s, " \n");
    			for (int i1 = 0; str.hasMoreTokens() && i1 < n; i1++) {
    				x[i1] = Double.parseDouble(str.nextToken());
    				y[i1] = Double.parseDouble(str.nextToken());
    			}
    					// 計算
    			int sw = least(m, n, x, y, z);
    
    			if (sw > 0)
    				ta2.setText("逆行列が存在しない");
    			else {
    				ta2.setText("");
    				for (int i1 = 0; i1 < m+1; i1++)
    					ta2.append((m-i1) + " 次の係数 " + String.format("%.5f",z[i1]) + "\n");
    			}
    		}
    	}
    
    	/*************************************/
    	/* 最初2乗法                        */
    	/*      m : 多項式の次数             */
    	/*      n : データの数               */
    	/*      x,y : データ                 */
    	/*      z : 多項式の係数(高次から) */
    	/*      return : =0 : 正常           */
    	/*               =1 : エラー         */
    	/*************************************/
    	int least(int m, int n, double x[], double y[], double z[])
    	{
    		double A[][], w[][], x1, x2;
    		int i1, i2, i3, sw = 0;
    
    		m++;
    		w = new double [m][m+1];
    		A = new double [n][m];
    
    		for (i1 = 0; i1 < n; i1++) {
    			A[i1][m-2] = x[i1];
    			A[i1][m-1] = 1.0;
    			x1         = A[i1][m-2];
    			x2         = x1;
    			for (i2 = m-3; i2 >= 0; i2--) {
    				x2        *= x1;
    				A[i1][i2]  = x2;
    			}
    		}
    
    		for (i1 = 0; i1 < m; i1++) {
    			for (i2 = 0; i2 < m; i2++) {
    				w[i1][i2] = 0.0;
    				for (i3 = 0; i3 < n; i3++)
    					w[i1][i2] += A[i3][i1] * A[i3][i2];
    			}
    		}
    
    		for (i1 = 0; i1 < m; i1++) {
    			w[i1][m] = 0.0;
    			for (i2 = 0; i2 < n; i2++)
    				w[i1][m] += A[i2][i1] * y[i2];
    		}
    
    		sw = gauss(w, m, 1, 1.0e-10);
    
    		if (sw == 0) {
    			for (i1 = 0; i1 < m; i1++)
    				z[i1] = w[i1][m];
    		}
    		else
    			sw = 1;
    
    		return sw;
    	}
    
    	/******************************************/
    	/* 線形連立方程式を解く(逆行列を求める) */
    	/*      w : 方程式の左辺及び右辺          */
    	/*      n : 方程式の数                    */
    	/*      m : 方程式の右辺の列の数          */
    	/*      eps : 逆行列の存在を判定する規準  */
    	/*      return : =0 : 正常                */
    	/*               =1 : 逆行列が存在しない  */
    	/******************************************/
    	int gauss(double w[][], int n, int m, double eps)
    	{
    		double y1, y2;
    		int ind = 0, nm, m1, m2, i1, i2, i3;
    
    		nm = n + m;
    
    		for (i1 = 0; i1 < n && ind == 0; i1++) {
    
    			y1 = .0;
    			m1 = i1 + 1;
    			m2 = 0;
    						// ピボット要素の選択
    			for (i2 = i1; i2 < n; i2++) {
    				y2 = Math.abs(w[i2][i1]);
    				if (y1 < y2) {
    					y1 = y2;
    					m2 = i2;
    				}
    			}
    						// 逆行列が存在しない
    			if (y1 < eps)
    				ind = 1;
    						// 逆行列が存在する
    			else {
    							// 行の入れ替え
    				for (i2 = i1; i2 < nm; i2++) {
    					y1        = w[i1][i2];
    					w[i1][i2] = w[m2][i2];
    					w[m2][i2] = y1;
    				}
    							// 掃き出し操作
    				y1 = 1.0 / w[i1][i1];
    
    				for (i2 = m1; i2 < nm; i2++)
    					w[i1][i2] *= y1;
    
    				for (i2 = 0; i2 < n; i2++) {
    					if (i2 != i1) {
    						for (i3 = m1; i3 < nm; i3++)
    							w[i2][i3] -= w[i2][i1] * w[i1][i3];
    					}
    				}
    			}
    		}
    
    		return(ind);
    	}
    }
    			

  3. JavaScript

      ここをクリックすると,任意のデータに対して画面上で実行することができます.
    <!DOCTYPE HTML>
    
    <HTML>
    
    <HEAD>
    
    	<TITLE>最小二乗法(多項式近似)</TITLE>
    	<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
    	<SCRIPT TYPE="text/javascript">
    		function main()
    		{
    					// データ
    			let m = parseInt(document.getElementById("m").value);
    			let n = parseInt(document.getElementById("n").value);
    			let x = new Array();
    			let y = new Array();
    			let z = new Array();
    			let s = (document.getElementById("data").value).split(/ {1,}|\n{1,}/);
    			let k = 0;
    			for (let i1 = 0; i1 < 2*n; i1 += 2) {
    				x[k] = parseFloat(s[i1]);
    				y[k] = parseFloat(s[i1+1]);
    				k++;
    			}
    					// 計算
    			let sw = least(m, n, x, y, z);
    
    			if (sw > 0)
    				document.getElementById("ans").value = "逆行列が存在しない";
    			else {
    				let str = "";
    				for (let i1 = 0; i1 < m+1; i1++)
    					str = str + (m-i1) + " 次の係数 " + z[i1] + "\n";
    				document.getElementById("ans").value = str;
    			}
    		}
    
    		/*************************************/
    		/* 最初2乗法                        */
    		/*      m : 多項式の次数             */
    		/*      n : データの数               */
    		/*      x,y : データ                 */
    		/*      z : 多項式の係数(高次から) */
    		/*      return : =0 : 正常           */
    		/*               =1 : エラー         */
    		/*************************************/
    		function least(m, n, x, y, z)
    		{
    			let x1;
    			let x2;
    			let i1;
    			let i2;
    			let i3;
    			let sw = 0;
    
    			m++;
    			let w = new Array();
    			for (i1 = 0; i1 < m; i1++)
    				w[i1] = new Array();
    			let A = new Array();
    			for (i1 = 0; i1 < n; i1++)
    				A[i1] = new Array();
    
    			for (i1 = 0; i1 < n; i1++) {
    				A[i1][m-2] = x[i1];
    				A[i1][m-1] = 1.0;
    				x1         = A[i1][m-2];
    				x2         = x1;
    				for (i2 = m-3; i2 >= 0; i2--) {
    					x2        *= x1;
    					A[i1][i2]  = x2;
    				}
    			}
    
    			for (i1 = 0; i1 < m; i1++) {
    				for (i2 = 0; i2 < m; i2++) {
    					w[i1][i2] = 0.0;
    					for (i3 = 0; i3 < n; i3++)
    						w[i1][i2] += A[i3][i1] * A[i3][i2];
    				}
    			}
    
    			for (i1 = 0; i1 < m; i1++) {
    				w[i1][m] = 0.0;
    				for (i2 = 0; i2 < n; i2++)
    					w[i1][m] += A[i2][i1] * y[i2];
    			}
    
    			sw = gauss(w, m, 1, 1.0e-10);
    
    			if (sw == 0) {
    				for (i1 = 0; i1 < m; i1++)
    					z[i1] = w[i1][m];
    			}
    			else
    				sw = 1;
    
    			return sw;
    		}
    
    		/******************************************/
    		/* 線形連立方程式を解く(逆行列を求める) */
    		/*      w : 方程式の左辺及び右辺          */
    		/*      n : 方程式の数                    */
    		/*      m : 方程式の右辺の列の数          */
    		/*      eps : 逆行列の存在を判定する規準  */
    		/*      return : =0 : 正常                */
    		/*               =1 : 逆行列が存在しない  */
    		/******************************************/
    		function gauss(w, n, m, eps)
    		{
    			let y1;
    			let y2;
    			let ind = 0;
    			let nm;
    			let m1;
    			let m2;
    			let i1;
    			let i2;
    			let i3;
    
    			nm = n + m;
    
    			for (i1 = 0; i1 < n && ind == 0; i1++) {
    
    				y1 = .0;
    				m1 = i1 + 1;
    				m2 = 0;
    						// ピボット要素の選択
    				for (i2 = i1; i2 < n; i2++) {
    					y2 = Math.abs(w[i2][i1]);
    					if (y1 < y2) {
    						y1 = y2;
    						m2 = i2;
    					}
    				}
    						// 逆行列が存在しない
    				if (y1 < eps)
    					ind = 1;
    						// 逆行列が存在する
    				else {
    							// 行の入れ替え
    					for (i2 = i1; i2 < nm; i2++) {
    						y1        = w[i1][i2];
    						w[i1][i2] = w[m2][i2];
    						w[m2][i2] = y1;
    					}
    							// 掃き出し操作
    					y1 = 1.0 / w[i1][i1];
    
    					for (i2 = m1; i2 < nm; i2++)
    						w[i1][i2] *= y1;
    
    					for (i2 = 0; i2 < n; i2++) {
    						if (i2 != i1) {
    							for (i3 = m1; i3 < nm; i3++)
    								w[i2][i3] -= w[i2][i1] * w[i1][i3];
    						}
    					}
    				}
    			}
    
    			return(ind);
    		}
    	</SCRIPT>
    
    </HEAD>
    
    <BODY STYLE="font-size: 130%; background-color: #eeffee;">
    
    	<H2 STYLE="text-align:center"><B>最小二乗法(多項式近似)</B></H2>
    
    	<DL>
    		<DT>  テキストフィールドおよびテキストエリアには,例として,テキストエリアに与えられた 10 個の点を直線で近似する場合に対する値が設定されています.他の問題を実行する場合は,それらを適切に修正してください.
    	</DL>
    
    	<DIV STYLE="text-align:center">
    		多項式の次数:<INPUT ID="m" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="1"> 
    		データの数:<INPUT ID="n" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="10"><BR><BR>
    		データ(x y):<TEXTAREA ID="data" COLS="30" ROWS="15" STYLE="font-size: 100%">
    0  2.450
    1  2.615
    2  3.276
    3  3.294
    4  3.778
    5  4.009
    6  3.920
    7  4.267
    8  4.805
    9  5.656
    		</TEXTAREA> 
    		<BUTTON STYLE="font-size: 100%; background-color: pink" onClick="main()">OK</BUTTON><BR><BR>
    		<TEXTAREA ID="ans" COLS="40" ROWS="5" STYLE="font-size: 100%;"></TEXTAREA>
    	</DIV>
    
    </BODY>
    
    </HTML>
    			

  4. PHP

    <?php
    
    /****************************/
    /* 最小二乗法(多項式近似) */
    /*      coded by Y.Suganuma */
    /****************************/
    
    	fscanf(STDIN, "%d %d", $m, $n);   // 多項式の次数とデータの数
    
    	$x = array($n);
    	$y = array($n);
    
    	for ($i1 = 0; $i1 < $n; $i1++)   // データ
    		fscanf(STDIN, "%lf %lf", $x[$i1], $y[$i1]);
    
    	$z = least($m, $n, $x, $y);
    
    	if ($z != NULL) {
    		printf("結果\n");
    		for ($i1 = 0; $i1 < $m+1; $i1++)
    			printf("   %d 次の係数 %f\n", $m-$i1, $z[$i1]);
    	}
    	else
    		printf("***error  逆行列を求めることができませんでした\n");
    
    /******************************************/
    /* 最初2乗法                             */
    /*      m : 多項式の次数                  */
    /*      n : データの数                    */
    /*      x,y : データ                      */
    /*      return : 多項式の係数(高次から) */
    /*               エラーの場合はNULLを返す */
    /******************************************/
    function least($m, $n, $x, $y)
    {
    	$m++;
    	$z = array($m);
    	$w = array($m);
    	for ($i1 = 0; $i1 < $m; $i1++)
    		$w[$i1] = array($m+1);
    	$A = array($n);
    
    	for ($i1 = 0; $i1 < $n; $i1++) {
    		$A[$i1]       = array($m);
    		$A[$i1][$m-2] = $x[$i1];
    		$A[$i1][$m-1] = 1.0;
    		$x1           = $A[$i1][$m-2];
    		$x2           = $x1;
    		for ($i2 = $m-3; $i2 >= 0; $i2--) {
    			$x2        *= $x1;
    			$A[$i1][$i2]  = $x2;
    		}
    	}
    
    	for ($i1 = 0; $i1 < $m; $i1++) {
    		for ($i2 = 0; $i2 < $m; $i2++) {
    			$w[$i1][$i2] = 0.0;
    			for ($i3 = 0; $i3 < $n; $i3++)
    				$w[$i1][$i2] += $A[$i3][$i1] * $A[$i3][$i2];
    		}
    	}
    
    	for ($i1 = 0; $i1 < $m; $i1++) {
    		$w[$i1][$m] = 0.0;
    		for ($i2 = 0; $i2 < $n; $i2++)
    			$w[$i1][$m] += $A[$i2][$i1] * $y[$i2];
    	}
    
    	$sw = gauss($w, $m, 1, 1.0e-10);
    
    	if ($sw == 0) {
    		for ($i1 = 0; $i1 < $m; $i1++)
    			$z[$i1] = $w[$i1][$m];
    	}
    	else
    		$z = NULL;
    
    	return $z;
    }
    
    /*******************************************************/
    /* 線形連立方程式を解く(逆行列を求める)              */
    /*      w : 方程式の左辺及び右辺                       */
    /*      n : 方程式の数                                 */
    /*      m : 方程式の右辺の列の数                       */
    /*      eps : 正則性を判定する規準                     */
    /*      return : =0 : 正常                             */
    /*               =1 : 逆行列が存在しない               */
    /*******************************************************/
    function gauss(&$w, $n, $m, $eps)
    {
    	$ind = 0;
    	$nm  = $n + $m;
    
    	for ($i1 = 0; $i1 < $n && $ind == 0; $i1++) {
    
    		$y1 = .0;
    		$m1 = $i1 + 1;
    		$m2 = 0;
    
    		for ($i2 = $i1; $i2 < $n; $i2++) {
    			$y2 = abs($w[$i2][$i1]);
    			if ($y1 < $y2) {
    				$y1 = $y2;
    				$m2 = $i2;
    			}
    		}
    
    		if ($y1 < $eps)
    			$ind = 1;
    
    		else {
    
    			for ($i2 = $i1; $i2 < $nm; $i2++) {
    				$y1          = $w[$i1][$i2];
    				$w[$i1][$i2] = $w[$m2][$i2];
    				$w[$m2][$i2] = $y1;
    			}
    
    			$y1 = 1.0 / $w[$i1][$i1];
    
    			for ($i2 = $m1; $i2 < $nm; $i2++)
    				$w[$i1][$i2] *= $y1;
    
    			for ($i2 = 0; $i2 < $n; $i2++) {
    				if ($i2 != $i1) {
    					for ($i3 = $m1; $i3 < $nm; $i3++)
    						$w[$i2][$i3] -= $w[$i2][$i1] * $w[$i1][$i3];
    				}
    			}
    		}
    	}
    
    	return($ind);
    }
    
    /*
    -----------データ例1(コメント部分を除いて下さい)---------
    1 10   // 多項式の次数とデータの数
    0  2.450   // 以下,(x, y)
    1  2.615
    2  3.276
    3  3.294
    4  3.778
    5  4.009
    6  3.920
    7  4.267
    8  4.805
    9  5.656
    
    -------------------------データ例2-------------------------
    2 6
    10  28.2
    15  47.0
    20  44.4
    26  32.8
    32  20.8
    40   0.8
    */
    
    ?>
    			

  5. Ruby

    ############################
    # 最小二乗法(多項式近似)
    #      coded by Y.Suganuma
    ############################
    
    ############################################
    # 線形連立方程式を解く(逆行列を求める)
    #      w : 方程式の左辺及び右辺
    #      n : 方程式の数
    #      m : 方程式の右辺の列の数
    #      eps : 逆行列の存在を判定する規準
    #      return : =0 : 正常
    #               =1 : 逆行列が存在しない
    #      coded by Y.Suganuma
    ############################################
    
    def gauss(w, n, m, eps)
    
    	nm  = n + m;
    	ind = 0
    
    	for i1 in 0 ... n
    
    		y1 = 0.0
    		m1 = i1 + 1
    		m2 = 0
    					# ピボット要素の選択
    		for i2 in i1 ... n
    			y2 = w[i2][i1].abs()
    			if y1 < y2
    				y1 = y2
    				m2 = i2
    			end
    		end
    					# 逆行列が存在しない
    		if y1 < eps
    			ind = 1
    			break
    					# 逆行列が存在する
    		else   # 行の入れ替え
    			for i2 in i1 ... nm
    				y1        = w[i1][i2]
    				w[i1][i2] = w[m2][i2]
    				w[m2][i2] = y1
    			end
    						# 掃き出し操作
    			y1 = 1.0 / w[i1][i1]
    
    			for i2 in m1 ... nm
    				w[i1][i2] *= y1
    			end
    
    			for i2 in 0 ... n
    				if i2 != i1
    					for i3 in m1 ... nm
    						w[i2][i3] -= (w[i2][i1] * w[i1][i3])
    					end
    				end
    			end
    		end
    	end
    
    	return ind
    end
    
    #########################################
    # 最小二乗法
    #      m : 多項式の次数
    #      n : データの数
    #      x,y : データ
    #      return : 多項式の係数(高次から)
    #               エラーの場合はNULLを返す
    #      coded by Y.Suganuma
    #########################################
    
    def least(m, n, x, y)
    
    	m += 1
    	z  = Array.new(m)
    	w  = Array.new(m)
    	for i1 in 0 ... m
    		w[i1] = Array.new(m+1)
    	end
    	a = Array.new(n)
    	for i1 in 0 ... n
    		a[i1] = Array.new(m)
    	end
    
    	for i1 in 0 ... n
    		a[i1][m-2] = x[i1]
    		a[i1][m-1] = 1.0
    		x1         = a[i1][m-2]
    		x2         = x1
    		i2         = m - 3
    		while i2 > -1
    			x2        *= x1
    			a[i1][i2]  = x2
    			i2        -= 1
    		end
    	end
    
    	for i1 in 0 ... m
    		for i2 in 0 ... m
    			w[i1][i2] = 0.0
    			for i3 in 0 ... n
    				w[i1][i2] += a[i3][i1] * a[i3][i2]
    			end
    		end
    	end
    
    	for i1 in 0 ... m
    		w[i1][m] = 0.0
    		for i2 in 0 ... n
    			w[i1][m] += a[i2][i1] * y[i2]
    		end
    	end
    
    	sw = gauss(w, m, 1, 1.0e-10)
    
    	if sw == 0
    		for i1 in 0 ... m
    			z[i1] = w[i1][m]
    		end
    	else
    		z = Array.new(0)
    	end
    
    	return z
    end
    
    s = gets().split(" ")
    m = Integer(s[0])   # 多項式の次数
    n = Integer(s[1])   # データの数
    
    x = Array.new(n)
    y = Array.new(n)
    
    for i1 in 0 ... n   # データ
    	s     = gets().split()
    	x[i1] = Float(s[0])
    	y[i1] = Float(s[1])
    end
    
    z = least(m, n, x, y)
    
    if z.size  > 0
    	print("結果\n")
    	for i1 in 0 ... m+1
    		print("   " + String(m-i1) + " 次の係数 " + String(z[i1]) + "\n")
    	end
    else
    	print("***error  逆行列を求めることができませんでした\n")
    end
    
    =begin
    ------------データ例1(コメント部分を除いて下さい)--------
    1 10   # 多項式の次数とデータの数
    0  2.450   # 以下,(x, y)
    1  2.615
    2  3.276
    3  3.294
    4  3.778
    5  4.009
    6  3.920
    7  4.267
    8  4.805
    9  5.656
    
    -------------------------データ例2-------------------------
    2 6
    10  28.2
    15  47.0
    20  44.4
    26  32.8
    32  20.8
    40   0.8
    =end
    			

  6. Python

    # -*- coding: UTF-8 -*-
    import numpy as np
    import sys
    from math import *
    
    ############################################
    # 線形連立方程式を解く(逆行列を求める)
    #      w : 方程式の左辺及び右辺
    #      n : 方程式の数
    #      m : 方程式の右辺の列の数
    #      eps : 逆行列の存在を判定する規準
    #      return : =0 : 正常
    #               =1 : 逆行列が存在しない
    #      coded by Y.Suganuma
    ############################################
    
    def gauss(w, n, m, eps) :
    
    	nm  = n + m
    	ind = 0
    
    	for i1 in range(0, n) :
    
    		y1 = 0.0
    		m1 = i1 + 1
    		m2 = 0
    					# ピボット要素の選択
    		for i2 in range(i1, n) :
    			y2 = abs(w[i2][i1])
    			if y1 < y2 :
    				y1 = y2
    				m2 = i2
    					# 逆行列が存在しない
    		if y1 < eps :
    			ind = 1
    			break
    					# 逆行列が存在する
    		else :   # 行の入れ替え
    			for i2 in range(i1, nm) :
    				y1        = w[i1][i2]
    				w[i1][i2] = w[m2][i2]
    				w[m2][i2] = y1
    						# 掃き出し操作
    			y1 = 1.0 / w[i1][i1]
    
    			for i2 in range(m1, nm) :
    				w[i1][i2] *= y1
    
    			for i2 in range(0, n) :
    				if i2 != i1 :
    					for i3 in range(m1, nm) :
    						w[i2][i3] -= (w[i2][i1] * w[i1][i3])
    
    	return ind
    
    #########################################
    # 最小二乗法
    #      m : 多項式の次数
    #      n : データの数
    #      x,y : データ
    #      return : 多項式の係数(高次から)
    #               エラーの場合はNULLを返す
    #      coded by Y.Suganuma
    #########################################
    
    def least(m, n, x, y) :
    
    	m += 1
    	z  = np.empty(m, np.float)
    	w  = np.empty((m, m+1), np.float)
    	A  = np.empty((n, m), np.float)
    
    	for i1 in range(0, n) :
    		A[i1][m-2] = x[i1]
    		A[i1][m-1] = 1.0
    		x1         = A[i1][m-2]
    		x2         = x1
    		for i2 in range(m-3, -1, -1) :
    			x2        *= x1
    			A[i1][i2]  = x2
    
    	for i1 in range(0, m) :
    		for i2 in range(0, m) :
    			w[i1][i2] = 0.0
    			for i3 in range(0, n) :
    				w[i1][i2] += A[i3][i1] * A[i3][i2]
    
    	for i1 in range(0, m) :
    		w[i1][m] = 0.0
    		for i2 in range(0, n) :
    			w[i1][m] += A[i2][i1] * y[i2]
    
    	sw = gauss(w, m, 1, 1.0e-10)
    
    	if sw == 0 :
    		for i1 in range(0, m) :
    			z[i1] = w[i1][m]
    	else :
    		z = np.empty(0, np.float)
    
    	return z
    
    ############################
    # 最小二乗法(多項式近似)
    #      coded by Y.Suganuma
    ############################
    
    line = sys.stdin.readline()
    s    = line.split()
    m    = int(s[0])   # 多項式の次数
    n    = int(s[1])   # データの数
    
    x    = np.empty(n, np.float)
    y    = np.empty(n, np.float)
    
    for i1 in range(0, n) :   # データ
    	line  = sys.stdin.readline()
    	s     = line.split()
    	x[i1] = float(s[0])
    	y[i1] = float(s[1])
    
    z = least(m, n, x, y)
    
    if z.size  > 0 :
    	print("結果")
    	for i1 in range(0, m+1) :
    		print("   " + str(m-i1) + " 次の係数 " + str(z[i1]))
    else :
    	print("***error  逆行列を求めることができませんでした")
    
    """
    ------------データ例1(コメント部分を除いて下さい)--------
    1 10   # 多項式の次数とデータの数
    0  2.450   # 以下,(x, y)
    1  2.615
    2  3.276
    3  3.294
    4  3.778
    5  4.009
    6  3.920
    7  4.267
    8  4.805
    9  5.656
    
    -------------------------データ例2-------------------------
    2 6
    10  28.2
    15  47.0
    20  44.4
    26  32.8
    32  20.8
    40   0.8
    """
    			

  7. C#

    /****************************/
    /* 最小二乗法(多項式近似) */
    /*      coded by Y.Suganuma */
    /****************************/
    using System;
    
    class Program
    {
    	static void Main()
    	{
    		Test1 ts = new Test1();
    	}
    }
    
    class Test1
    {
    	public Test1()
    	{
    		char[] charSep = new char[] {' '};
    					// 多項式の次数とデータの数
    		String[] str = Console.ReadLine().Split(charSep, StringSplitOptions.RemoveEmptyEntries);
    		int m        = int.Parse(str[0]);
    		int n        = int.Parse(str[1]);
    
    		double[] x = new double [n];
    		double[] y = new double [n];
    		double[] z = new double [m+1];
    					// データ
    		for (int i1 = 0; i1 < n; i1++) {
    			str   = Console.ReadLine().Split(charSep, StringSplitOptions.RemoveEmptyEntries);
    			x[i1] = double.Parse(str[0]);
    			y[i1] = double.Parse(str[1]);
    		}
    
    		int sw = least(m, n, x, y, z);
    
    		if (sw == 0) {
    			Console.WriteLine("結果");
    			for (int i1 = 0; i1 < m+1; i1++)
    				Console.WriteLine("   " + (m-i1) + " 次の係数 " + z[i1]);
    		}
    		else
    			Console.WriteLine("***error  逆行列を求めることができませんでした");
    	}
    
    	/*************************************/
    	/* 最初2乗法                        */
    	/*      m : 多項式の次数             */
    	/*      n : データの数               */
    	/*      x,y : データ                 */
    	/*      z : 多項式の係数(高次から) */
    	/*      return : =0 : 正常           */
    	/*               =1 : エラー         */
    	/*************************************/
    	int least(int m, int n, double[] x, double[] y, double[] z)
    	{
    		m++;
    		double[][] w = new double [m][];
    		for (int i1 = 0; i1 < m; i1++)
    			w[i1] = new double[m+1];
    		double[][] A = new double [n][];
    		for (int i1 = 0; i1 < n; i1++)
    			A[i1] = new double[m];
    
    		for (int i1 = 0; i1 < n; i1++) {
    			A[i1][m-2] = x[i1];
    			A[i1][m-1] = 1.0;
    			double x1  = A[i1][m-2];
    			double x2  = x1;
    			for (int i2 = m-3; i2 >= 0; i2--) {
    				x2        *= x1;
    				A[i1][i2]  = x2;
    			}
    		}
    
    		for (int i1 = 0; i1 < m; i1++) {
    			for (int i2 = 0; i2 < m; i2++) {
    				w[i1][i2] = 0.0;
    				for (int i3 = 0; i3 < n; i3++)
    					w[i1][i2] += A[i3][i1] * A[i3][i2];
    			}
    		}
    
    		for (int i1 = 0; i1 < m; i1++) {
    			w[i1][m] = 0.0;
    			for (int i2 = 0; i2 < n; i2++)
    				w[i1][m] += A[i2][i1] * y[i2];
    		}
    
    		int sw = gauss(w, m, 1, 1.0e-10);
    
    		if (sw == 0) {
    			for (int i1 = 0; i1 < m; i1++)
    				z[i1] = w[i1][m];
    		}
    		else
    			sw = 1;
    
    		return sw;
    	}
    
    	/*******************************************************/
    	/* 線形連立方程式を解く(逆行列を求める)              */
    	/*      w : 方程式の左辺及び右辺                       */
    	/*      n : 方程式の数                                 */
    	/*      m : 方程式の右辺の列の数                       */
    	/*      eps : 逆行列の存在を判定する規準               */
    	/*      return : =0 : 正常                             */
    	/*               =1 : 逆行列が存在しない               */
    	/*******************************************************/
    	int gauss(double[][] w, int n, int m, double eps)
    	{
    		int ind = 0;
    		int nm  = n + m;
    
    		for (int i1 = 0; i1 < n && ind == 0; i1++) {
    
    			double y1 = .0;
    			int m1    = i1 + 1;
    			int m2    = 0;
    						// ピボット要素の選択
    			for (int i2 = i1; i2 < n; i2++) {
    				double y2 = Math.Abs(w[i2][i1]);
    				if (y1 < y2) {
    					y1 = y2;
    					m2 = i2;
    				}
    			}
    						// 逆行列が存在しない
    			if (y1 < eps)
    				ind = 1;
    						// 逆行列が存在する
    			else {
    							// 行の入れ替え
    				for (int i2 = i1; i2 < nm; i2++) {
    					y1        = w[i1][i2];
    					w[i1][i2] = w[m2][i2];
    					w[m2][i2] = y1;
    				}
    							// 掃き出し操作
    				y1 = 1.0 / w[i1][i1];
    
    				for (int i2 = m1; i2 < nm; i2++)
    					w[i1][i2] *= y1;
    
    				for (int i2 = 0; i2 < n; i2++) {
    					if (i2 != i1) {
    						for (int i3 = m1; i3 < nm; i3++)
    							w[i2][i3] -= w[i2][i1] * w[i1][i3];
    					}
    				}
    			}
    		}
    
    		return ind;
    	}
    }
    
    /*
    -----------データ例1(コメント部分を除いて下さい)---------
    1 10   // 多項式の次数とデータの数
    0  2.450   // 以下,(x, y)
    1  2.615
    2  3.276
    3  3.294
    4  3.778
    5  4.009
    6  3.920
    7  4.267
    8  4.805
    9  5.656
    
    -------------------------データ例2-------------------------
    2 6
    10  28.2
    15  47.0
    20  44.4
    26  32.8
    32  20.8
    40   0.8
    */
    			

  8. VB

    '**************************'
    ' 最小二乗法(多項式近似) '
    '      coded by Y.Suganuma '
    '**************************'
    Imports System.Text.RegularExpressions
    
    Module Test
    
    	Sub Main()
    
    		Dim MS As Regex = New Regex("\s+") 
    					' 多項式の次数とデータの数
    		Dim str() As String = MS.Split(Console.ReadLine().Trim())
    		Dim m As Integer    = Integer.Parse(str(0))
    		Dim n As Integer    = Integer.Parse(str(1))
    
    		Dim x(n) As Double
    		Dim y(n) As Double
    		Dim z(m+1) As Double
    					' データ
    		For i1 As Integer = 0 To n-1
    			str   = MS.Split(Console.ReadLine().Trim())
    			x(i1) = Double.Parse(str(0))
    			y(i1) = Double.Parse(str(1))
    		Next
    
    		Dim sw As Integer = least(m, n, x, y, z)
    
    		If sw = 0
    			Console.WriteLine("結果")
    			For i1 As Integer = 0 To m
    				Console.WriteLine("   " & (m-i1) & " 次の係数 " & z(i1))
    			Next
    		Else
    			Console.WriteLine("***error  逆行列を求めることができませんでした")
    		End If
    
    	End Sub
    
    	'***********************************'
    	' 最初2乗法                        '
    	'      m : 多項式の次数             '
    	'      n : データの数               '
    	'      x,y : データ                 '
    	'      z : 多項式の係数(高次から) '
    	'      return : =0 : 正常           '
    	'               =1 : エラー         '
    	'***********************************'
    	Function least(m As Integer, n As Integer, x() As Double, y() As Double, z() As Double)
    
    		m += 1
    		Dim w(m,m+1) As Double
    		Dim A(n,m) As Double
    
    		For i1 As Integer = 0 To n-1
    			A(i1,m-2) = x(i1)
    			A(i1,m-1) = 1.0
    			Dim x1 As Double  = A(i1,m-2)
    			Dim x2 As Double  = x1
    			For i2 As Integer = m-3 To 0 Step -1
    				x2       *= x1
    				A(i1,i2)  = x2
    			Next
    		Next
    
    		For i1 As Integer = 0 To m-1
    			For i2 As Integer = 0 To m-1
    				w(i1,i2) = 0.0
    				For i3 As Integer = 0 To n-1
    					w(i1,i2) += A(i3,i1) * A(i3,i2)
    				Next
    			Next
    		Next
    
    		For i1 As Integer = 0 To m-1
    			w(i1,m) = 0.0
    			For i2 As Integer = 0 To n-1
    				w(i1,m) += A(i2,i1) * y(i2)
    			Next
    		Next
    
    		Dim sw As Integer = gauss(w, m, 1, 1.0e-10)
    
    		If sw = 0
    			For i1 As Integer = 0 To m-1
    				z(i1) = w(i1,m)
    			Next
    		Else
    			sw = 1
    		End If
    
    		Return sw
    
    	End Function
    
    	''''''''''''''''''''''''''''''''''''''''''
    	' 線形連立方程式を解く(逆行列を求める) '
    	'      w : 方程式の左辺及び右辺          '
    	'      n : 方程式の数                    '
    	'      m : 方程式の右辺の列の数          '
    	'      eps : 逆行列の存在を判定する規準  '
    	'      return : =0 : 正常                '
    	'               =1 : 逆行列が存在しない  '
    	''''''''''''''''''''''''''''''''''''''''''
    	Function gauss(w(,) As Double, n As Integer, m As Integer, eps As Double) As Integer
    
    		Dim ind As Integer = 0
    		Dim nm As Integer  = n + m
    
    		Dim i1 As Integer = 0
    		Do While i1 < n and ind = 0
    
    			Dim y1 As Double  = 0.0
    			Dim m1 As Integer = i1 + 1
    			Dim m2 As Integer = 0
    						' ピボット要素の選択
    			For i2 As Integer = i1 To n-1
    				Dim y2 As Double = Math.Abs(w(i2,i1))
    				If y1 < y2
    					y1 = y2
    					m2 = i2
    				End If
    			Next
    						' 逆行列が存在しない
    			If y1 < eps
    				ind = 1
    						' 逆行列が存在する
    			Else
    							' 行の入れ替え
    				For i2 As Integer = i1 To nm-1
    					y1       = w(i1,i2)
    					w(i1,i2) = w(m2,i2)
    					w(m2,i2) = y1
    				Next
    							' 掃き出し操作
    				y1 = 1.0 / w(i1,i1)
    
    				For i2 As Integer = m1 To nm-1
    					w(i1,i2) *= y1
    				Next
    
    				For i2 As Integer = 0 To n-1
    					If i2 <> i1
    						For i3 As Integer = m1 To nm-1
    							w(i2,i3) -= w(i2,i1) * w(i1,i3)
    						Next
    					End If
    				Next
    			End If
    			i1 += 1
    		Loop
    
    		Return ind
    
    	End Function
    
    End Module
    
    /*
    -----------データ例1(コメント部分を除いて下さい)---------
    1 10   // 多項式の次数とデータの数
    0  2.450   // 以下,(x, y)
    1  2.615
    2  3.276
    3  3.294
    4  3.778
    5  4.009
    6  3.920
    7  4.267
    8  4.805
    9  5.656
    
    -------------------------データ例2-------------------------
    2 6
    10  28.2
    15  47.0
    20  44.4
    26  32.8
    32  20.8
    40   0.8
    */
    			

静岡理工科大学 菅沼ホーム 目次 索引