a*(1-exp(-bx))+cの近似の方法

質問

a*(1-exp(-bx))+cの近似の方法

x　　y
0　　1.00
2　　1.90
5　　2.96
7　　3.51
10　　4.16
20　　5.32
30　　5.75

ある実験で上記のようなデータを取得しました。
ほとんどのyデータはx=30で3～30の間でほぼ飽和しています。
cは0～5の間を取ります。

現在市販のグラフソフトを使用してa*(1-exp(-bx))+cで近似曲線を得ています。
またデータを取るソフトはVCで作成しています。

データ数を毎回、グラフソフトを立ち上げてひとつづつ近似するのは
手間がかかってしまいますし、データ数も膨大になってきました。

そこで作成しているソフト内に組み込もうと思うのですが、どういった方法が適切なのかが
検討がつきません。方法またはそのようなことのできるライブラリなどご存知であればご回答よろしくお願いします。

回答No.3

和泉 博（@hiroshi09s） · ベストアンサー

GnuPlotは独立したひとつのプログラムですから、定数a,b,c を直接利用することはできません。fit.logファイルを読むしかありません。
　ソースは2000字制限で全部を載せることができないため、#2の「 input_str();if () else 」文をコピペしてください。

/* Gcc on Mac OSX
* file name: r_aso.c
* compile: gcc r_aso.c
* execution: ./a.out
*/

#include <stdio.h>
#include <stdlib.h> //exit(),system()
#include <string.h> //strtok(),strstr()

#define SIZE 128
#define FILE_NAME "jikken.dat"
#define GP_LOGFILE "fit.log"
#define SPC " "
#define KEYWORD "+/-"
#define input_str(c,x) {printf(c);fgets(x,SIZE,stdin);}
#define set_token(v,x) {strtok(x,SPC);strtok(NULL,SPC);v=atof(strtok(NULL,SPC));}

int main(void) {
int i;
char buff[SIZE];
FILE *fp,*gp;
double a,b,c;

/*ここに
input_str("Data input? (y/n) ",buff);
if (*buff=='y' || *buff=='Y') {
....
をコピペのこと。 */

// ***** GnuPlot *****
if ((gp=popen("gnuplot","w"))==NULL) {
perror("popen");
exit(1);
}
// 既存fit.logファイルの削除
sprintf(buff,"rm %s\n",GP_LOGFILE);
system(buff);
// fit.log の作成
fprintf(gp,"f(x)=a*(1-exp(-b*x))+c \n");
fprintf(gp,"a=5;b=0.1;c=1 \n");
fprintf(gp,"fit f(x) '%s' via a,b,c \n",FILE_NAME);
pclose(gp); //←ここでいったん GnuPlotを閉じる

//定数a,b,c の読み取り
if ((fp=fopen(GP_LOGFILE,"r"))==NULL) {
perror("fopen(GP_LOGFILE)");
exit(1);
}
fgets(buff,SIZE,fp);
while (strstr(buff,KEYWORD)==NULL){
fgets(buff,SIZE,fp);
}
set_token(a,buff);
fgets(buff,SIZE,fp);
set_token(b,buff);
fgets(buff,SIZE,fp);
set_token(c,buff);
fclose(fp);
sprintf(buff,"%f*(1-exp(-%f*x))+%f",a,b,c);

//GnuPlotによる描画
gp=popen("gnuplot","w");
fprintf(gp,"set xrange [0:40] \n");
fprintf(gp,"set yrange [0:8] \n");
fprintf(gp,"set title 'An example' \n");
fprintf(gp,"set xlabel 'X' \n");
fprintf(gp,"set ylabel 'f(x)' \n");
fprintf(gp,"plot '%s' with points pointtype 6 pointsize 1.5 \n",FILE_NAME);
fprintf(gp,"f(x)=%s \n", buff);
fprintf(gp,"replot f(x) with lines title 'f(x)= %s' \n", buff);
pclose(gp);

//結果のターミナル出力
printf("\n");
printf("GnuPlot fitting result:\n");
printf("\tf(x) = %s\n", buff);

return 0;
}

お礼

どうにか最小プログラムで実現できました。

ソフトに組み込んでみます。
本当にありがとうございました。

処理が重くなったときはそのときにほかの方法を考えます。

gnuplotって便利ですね。

gnuplotについてはこれから勉強していきます。

回答No.2

和泉 博（@hiroshi09s）

当方、Windows のVCではないMac OSX ですが、参考になれば幸いです。
プログラムは、日々データを加えて図化してはその傾向を視覚確認することを簡略化するのに重宝します。Windows版GnuPlotはMac OSXと使い方に細かな違いがあると思われますので「gnuplot>」のコマンドラインからfprintf(gp, ...) の内容を実行して確認してみてください。

/* Gcc on Mac OSX
* See http://t16web.lanl.gov/Kawano/gnuplot/misc2.html#9.3
* file name: raso.c
* compile: gcc raso.c
* execution: ./a.out
*/

#include <stdio.h>
#include <stdlib.h> //exit()
#define FILE_NAME "jikken.dat"
#define GP_LOGFILE "fit.log"
#define SIZE 128
#define input_str(c,x) {printf(c); fgets(x, SIZE, stdin);}

int main(void) {
int i;
char buff[SIZE];
FILE *fp, *gp;

input_str("Data input? (y/n) ", buff);
if (*buff == 'y' || *buff == 'Y') {
if ((fp = fopen(FILE_NAME, "a+")) == NULL) {
perror("fopen");
exit(1);
}

printf("\tData inputting and curve fitting.\n");
printf("\tLoop end: Push [return] key only.\n");
input_str("Input comment: ", buff);
fprintf(fp, "\n");
fprintf(fp, "# %s", buff);

i = 1;
printf("%d_ ", i);
input_str("x y: ", buff);
while (*buff != '\n') {
fprintf(fp, "%s", buff);
i++;
printf("%d_ ", i);
input_str("x y: ", buff);
}
fclose(fp);

} else {
if (!(*buff == 'n' || *buff == 'N')) {
printf("select y/n error.\n");
exit(1);
}
if ((fp = fopen(FILE_NAME, "r")) == NULL) {
perror("fopen");
exit(1);
}
printf("\tCurve fitting only.\n");
fclose(fp);
}

// GnuPlot
if ((gp = popen("gnuplot", "w")) == NULL) {
perror("popen");
exit(1);
}
// 既存fit.logファイルの削除
sprintf(buff, "rm %s\n", GP_LOGFILE);
system(buff);
// GnuPlot の実行
fprintf(gp, "f(x)=a*(1-exp(-b*x))+c \n");
fprintf(gp, "a=5; b=0.5; c=1 \n");
fprintf(gp, "fit f(x) '%s' via a,b,c \n", FILE_NAME);
fprintf(gp, "set title 'An example' \n");
fprintf(gp, "set xlabel 'X' \n");
fprintf(gp, "set ylabel 'Y' \n");
fprintf(gp, "plot f(x) with lines title 'a*(1-exp(-bx))+c', '%s' with points \n", FILE_NAME);
pclose(gp);
// 新fit.logファイルの閲覧
sprintf(buff, "more %s\n", GP_LOGFILE);
system(buff);

return 0;
}

お礼

gnuplotはきいたことがありましたが使ったことがなかったので、
Cから使えるって発想がなかったです。

なんか簡単にできそうな気がしてきました。試してみます。

補足

gnuplotから得られたパラメータ"a","b","c"をCの変数に格納できるでしょうか。
gnuplotの"print a,b,c"の値が取得できれば解決しそうです。

fit.logから探すとなると処理が多くなるとなんとも気持ちがわるいので、
どうにかならないかなと思っています。

もしお分かりであれば、ご回答よろしくお願いいたします。

回答No.1

kmee

近似式の形が決っているのなら、最小二乗法を使うのが一般的でしょう。
よく使われるものなので、解法の解説も多く、自分でプログラムを組んでもいいのです。
が、数値演算ライブラリに入っていることも多いので、それを使うと便利でしょう。

VCの例があるということで、こんなのはどうでしょうか
http://www.mlab.ice.uec.ac.jp/~ej-sib/numerical/numerical_levmar.html

お礼

こんな便利なものがあるものですね。

勉強してみます。
ご回答ありがとうございました。