最小二乗法

フォーラム(掲示板)ルール
フォーラム(掲示板)ルールはこちら  ※コードを貼り付ける場合は [code][/code] で囲って下さい。詳しくはこちら
たつみ

最小二乗法

#1

投稿記事 by たつみ » 12年前

何時間も考えてわけわからなくなりました。
手助けお願いします。
問題は x=(-4,-3,-2,-1,0,1,2,3,4,5,6)
y=(13.4,7.2,2.7,1.2,0.6,3.1,6.9,12.5,20.8,31.7,44.0)
を観測点とし、それに最もよく近似すると思われるf(x)を求めることです。
ガウス・ヨルダン法を使って解きます。
ガウス・ヨルダンのプログラムは確実にあっています。
わからないのはそれに渡す配列の作り方。
main関数。
結果をエクセルのグラフで表すこと。
作ったエラーだらけのソース貼ります。
#include<stdio.h>
#include<math.h>
#include<stdlib.h>

#define ERROR -1	//エラーリターン値
#define NORMAL 0	//ノーマルリターン値
#define mabs(x)  ( ( (x) > 0 ) ? (x):-(x))	//絶対値を返すマクロ
#define DBL_EPSILON 0.000000001


//Gauss_Jに渡す配列を作る↓

double setmatrixa(double x,int m,int n){
	
	int i,j,k;
	double a[100];
	
	for(i=0;i<m;i++){
		
		for(j=0;j<m;j++){
			
			for(k=0;k<n;k++){
				
				a[(m+2)*i+j]+=Mypow(x[k],(i+j));
				
				}
			}
		}
	
	return(a);
	
	}

//Gauss_Jに渡す配列を作る↓

double setmatrixb(double y,int m,int n){
		
		int i,k;
		double b[100];
		
		for(i=0;i<m;i++){
			
			for(k=0;k<n;k++){
				
				b+=y[k]*Mypow(x[k],i);
				
				}
			
			}
		
	return(b);
		
		}

//POW関数だとなぜかエラーが出るため、代用の関数↓

double Mypow(double x,int s){
			
			double buf_x;
			int i;
			buf_x=x;
			
			if(s==0){
				
				x=1;
				
				}
			else{
			
			for(i=0;i<s-1;i++){
				
				x=x*buf_x;
				
				}
			
			}
	return(x);
			
			}

//ここから


int Gauss_J(int n,double *a,double *b);
int partial(int n,int p,double *a,double *b);
void mswap(double *a,double *b);

int Gauss_J(int n,double *a,double *b){
	
	int p,i,j,l;	//カウンタ
	double pivot,c;	//ピボット値
	
	for(p=0;p<n;p++){	//1行目からn行目まで繰り返す
		
		if((partial(n,p,(double *)a,(double *)b))==ERROR){
			
			return(ERROR);
			
			}
		pivot=a[p*n+p];	//ピボットを取得
		for(i=p;i<n;i++){	//p行目のp列からn列まで
			
			a[p*n+i]/=pivot;	//係数行列のp行をpivotで割る
			
			}
		b[p]/=pivot;	//定数ベクトルのp行をpivotで割る
		
		for(l=0;l<n;l++){	//1行目からn行目まで
			
			if(l!=p){	//p行を除いて
				
				c=a[l*n+p];	//掃きだす
				for(j=p;j<n;j++){
					
					a[l*n+j]-=c*a[p*n+j];
					
					}
				b[[/url]-=c*b[p];
				
				}
			
			}
		}
	return(NORMAL);
	
	}

int partial(int n,int p,double *a,double *b){
		
		int l,m,max;	//l:カウンタ max:絶対値が最大値の係数を持つ行
		double pivot;	//ピボット
		
		for(l=max=p,pivot=a[p*n+p];l<n;l++){	//n行目まででp列の絶対値が最大のものを選ぶ
			
			if(mabs(pivot)<mabs(a[l*n+p])){
				
				max=l;
				pivot=a[l*n+p];
				
				}
			
			}
		if(mabs(pivot)<DBL_EPSILON){	//選ばれたピボットが0ならERROR(-1)を返す
			
			return(ERROR);
			
			}
		if(p!=max){	//選ばれた行がp行でなければ、行の交換を行う
			
			for(m=0;m<n;m++){	//係数行列のpとmaxを交換
				
				mswap(&a[p*n+m],&a[max*n+m]);
				
				}
			mswap(&b[p],&b[max]);	//定数のp行とmax行を交換
			
			}
		return (NORMAL);
		}

void mswap(double *a,double *b){
			
	double temp;
	
	temp=*a;
	*a=*b;
	*b=temp;
	
	}


//ここまで教科書丸写し



void main(void){
	
	double x[100],y[100],a[100],b[100];
	int i,n,m;
	int r;
	
	printf("データ数 = ");
	scanf("%d",&n);
	printf("%d\n",n);
	
	printf("次元数 = ");
	scanf("%d",&m);
	printf("%d\n",m);
	
	for(i=0;i<n;i++){
		
		printf("x = ");
		scanf("%lf",&x);
		printf("%f",x);
		
		printf("y = ");
		scanf("%lf",&y);
		printf("%f",y);
		
		}
		
		a=setmatrixa(x,m,n);
		b=setmatrixb(y,m,n);
		
		r=Gauss_J(m,a,b);
		
	if(r==0){	//エラー判定
		
		printf("isw = %d\n",r);
		
		for(i=0;i<m;i++){
		
		printf("X%d = %lf\n",i+1,b);	//結果出力
		
		}
		
		}
	
	else{
			
			printf("isw = %d",r);
			
			}
	}

box

Re:最小二乗法

#2

投稿記事 by box » 12年前

> 問題は x=(-4,-3,-2,-1,0,1,2,3,4,5,6)
> y=(13.4,7.2,2.7,1.2,0.6,3.1,6.9,12.5,20.8,31.7,44.0)
> を観測点とし、それに最もよく近似すると思われるf(x)を求めることです。

f(x)の式を求めたいのですか?
f(x)に条件はないのですか?例えば2次式であるとか。

あ、それとも、(x,y)の組が11個あるから、
xの10次式の係数(11個)を、11元連立一次方程式を立てて求めなさい、
ということなのでしょうか?

box

Re:最小二乗法

#3

投稿記事 by box » 12年前

> あ、それとも、(x,y)の組が11個あるから、
> xの10次式の係数(11個)を、11元連立一次方程式を立てて求めなさい、
> ということなのでしょうか?

与えられている条件だけでは、この考え方は無理っぽいですか…。
よくわからないですけど。

管理人

Re:最小二乗法

#4

投稿記事 by 管理人 » 12年前

>わからないのはそれに渡す配列の作り方。

計算内容は拝見しておりませんが、渡す配列の渡し方を知りたいならそこだけ個別に聞いたほうが聞きたい内容の情報が手に入るのではないでしょうか。

行いたいことの実装の仕方を質問されてはいかがでしょうか。

初心者

おしえてください

#5

投稿記事 by 初心者 » 12年前

はじめまして。大学に通っている者なのですが次のような問題がでてわからずに困っています。

「問題」
ポインタ、構造体、関数分割、ファイル分割のすべてを使ったプログラムをつくりなさい

という問題がでました。僕はIF文かswitch関数を使えば簡単なプログラムができるのではないかと思うのですがどのようにプログラムを組み立てればよいのかわかりません。もしよろしければ教えていただけないでしょうか?

バグ

Re:おしえてください

#6

投稿記事 by バグ » 12年前

過去レスで全く同じ質問をされていた方が居ましたが、その方は関数分割とファイル分割が分からないという事で質問されていました。あなたは何が分かりませんか?
せめて自分が何を作りたいかくらいは、決めてもらわないと誰もアドバイスのしようがないですよ。

初心者

Re:おしえてください

#7

投稿記事 by 初心者 » 12年前

すいません。授業を体調不良などで休んでいて友達に授業内容を聞いたりしたのですが今回の問題のことについてはあまり理解できていません。ポインタ、構造体、関数分割、ファイル分割のすべてを使ったプログラムで計算機のようなものは作れないでしょうか?

閉鎖

“C言語何でも質問掲示板” へ戻る