| 情報学部 | 菅沼ホーム | 目次 | 索引 |
/****************************/
/* 最小二乗法（多項式近似） */
/* 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;
}
/******************************************/
/* 最初２乗法 */
/* 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 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 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.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 逆行列を求めることができませんでした");
}
/*************************************/
/* 最初２乗法 */
/* 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 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 6
10 28.2
15 47.0
20 44.4
26 32.8
32 20.8
40 0.8
*/
<!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;
}
}
/*************************************/
/* 最初２乗法 */
/* 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>
<?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");
/******************************************/
/* 最初２乗法 */
/* 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 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 6
10 28.2
15 47.0
20 44.4
26 32.8
32 20.8
40 0.8
*/
?>
############################
# 最小二乗法（多項式近似）
# 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 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 6
10 28.2
15 47.0
20 44.4
26 32.8
32 20.8
40 0.8
=end
# -*- 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 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 6
10 28.2
15 47.0
20 44.4
26 32.8
32 20.8
40 0.8
"""
/****************************/
/* 最小二乗法（多項式近似） */
/* 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 逆行列を求めることができませんでした");
}
/*************************************/
/* 最初２乗法 */
/* 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 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 6
10 28.2
15 47.0
20 44.4
26 32.8
32 20.8
40 0.8
*/
'**************************'
' 最小二乗法（多項式近似） '
' 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
'***********************************'
' 最初２乗法 '
' 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 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 6
10 28.2
15 47.0
20 44.4
26 32.8
32 20.8
40 0.8
*/
| 情報学部 | 菅沼ホーム | 目次 | 索引 |