|
楼主 |
发表于 2010-6-11 15:39:27
|
显示全部楼层
找到了例子中的原型,发现其实很简单,但是初值的选择还是有一定的难度,网上资料很少,1stOpt具有全局最优的能力,但是看不到,希望哪位高手能提供一个初值的粗略估算方法
以下是实现的几个简单曲线拟合实例,难点就在于初值选择
#include <stdio.h>
#include <process.h>
//#include <stdlib.h>
#include <math.h>
//#include <float.h>
#include "levmar.h"
#include "compiler.h"
#include "misc.h"
#include "test1.h"
//因为程序是C++,而CLAPACK是f2c程序转换的C语言版本,所以在此处用extern关键字调用
extern"C"
{
#include <f2c.h>
#include <clapack.h>
}
#define ROSD 105.0
#define EV 2.71828183
double ava[6]={1,2,3,4,5,6};
double bva[6]={1,2,5,7,40,43};
///* Rosenbrock function, global minimum at (1, 1) */
//void ros(double *p, double *x, int m, int n, void *data)
//{
// register int i;
//
// for(i=0; i<n; ++i)
// x=((1.0-p[0])*(1.0-p[0]) + ROSD*(p[1]-p[0]*p[0])*(p[1]-p[0]*p[0]));
//}
//
//void jacros(double *p, double *jac, int m, int n, void *data)
//{
// register int i, j;
//
// for(i=j=0; i<n; ++i){
// jac[j++]=(-2 + 2*p[0]-4*ROSD*(p[1]-p[0]*p[0])*p[0]);
// jac[j++]=(2*ROSD*(p[1]-p[0]*p[0]));
// }
//}
//
//#define SIZE 4
//
//void func(double *p, double *x, int m, int n, void *data)
//{
// int i;
// for(i = 0; i < n; ++i)
// x = (p[0] - 3.1415926) * (p[0] - 3.1415926) +
// (p[1] - 2.0) * (p[1] - 2.0);
//}
/* Wood's function, minimum at (1, 1, 1, 1) */
//void wood(double *p, double *x, int m, int n, void *data)
//{
// register int i;
//
// for(i=0; i<n; i+=6)
// {
// x=10.0*(p[1] - p[0]*p[0]);
// x[i+1]=1.0 - p[0];
// x[i+2]=sqrt(90.0)*(p[3] - p[2]*p[2]);
// x[i+3]=1.0 - p[2];
// x[i+4]=sqrt(10.0)*(p[1]+p[3] - 2.0);
// x[i+5]=(p[1] - p[3])/sqrt(10.0);
// }
//}
//void func(double *p, double *x, int m, int n, void *data)
//{
// int i;
// double a[3]={1,2,3};
// double b[3];
// for(i=0;i<3;i++)
// b=2*a*a*a+3*a*a+a;
// for(i = 0; i < n; i+=3)
// {
// x = b-(p[0]*a*a*a+p[1]*a*a+p[2]*a);
// x[i+1] =b[i+1]-(p[0]*a[i+1]*a[i+1]*a[i+1]+p[1]*a[i+1]*a[i+1]+p[2]*a[i+1]);
// x[i+2] =b[i+2]-(p[0]*a[i+2]*a[i+2]*a[i+2]+p[1]*a[i+2]*a[i+2]+p[2]*a[i+2]);
// }
//}
//void jacfunc(double *p, double *jac, int m, int n, void *data)
//{
// register int i, j;
//
// for(i=j=0; i<n; ++i){
// jac[j++]=-1;
// jac[j++]=-1;
//
// jac[j++]=-4;
// jac[j++]=-2;
// }
//}
//void funcs(float x, float a[], float *y, float dyda[], int na)
//{
// int i;
// float fac,ex,arg;
// *y=0.0;
// for(i=1;i<=na-1;i+=3)
// {
// fac=a+a[i+1]*x;
// ex=exp(fac);
// arg=(1+ex)*(1+ex);
// *y=a[i+2]/(1+ex)+a[4];
// dyda=(-a[i+2])*ex/arg;
// dyda[i+1]=(-a[i+2])*x*ex/arg;
// dyda[i+2]=1/(1+ex);
// dyda[i+3]=1;
// }
//}
void Weibull(double *p, double *x, int m, int n, void *data)
{
int i;
//for(i=0;i<4;i++)
//bva=1-3*exp(-4*pow(ava,1));
//for(i=0; i<m; ++i) printf("%.7g ", b);printf("\n");
for(i = 0; i < n; i++)
{
double u = pow(ava,p[3]);
x = bva-(p[0]-p[1]*exp(-p[2]*u));
}
}
////void jacWeibull(double *p, double *jac, int m, int n, void *data)
////{
//// register int i, j;
////
//// for(i=j=0; i<n; ++i)
//// {
//// jac[j++]=-1;
//// jac[j++]=exp(-p[2]*pow(ava,p[3]));
//// jac[j++]=-p[1]*exp(-p[2]*pow(ava,p[3]))*-p[2]*pow(ava,p[3])*(log(ava)/log(EV));
//// }
////}
void Logistic(double *p, double *x, int m, int n, void *data)
{
int i;
//for(i=0;i<3;i++)
//b=2/(1+1*exp(-2*a));
//for(i=0; i<m; ++i) printf("%.7g ", b);printf("\n");
for(i = 0; i < n; i++)
{
x = bva-p[0]/(1+p[1]*exp(-p[2]*ava));
}
}
void funa(double *p, double *x, int m, int n, void *data)
{
int i;
//for(i=0; i<m; ++i) printf("%.7g ", b);printf("\n");
for(i = 0; i < n; i++)
{
x = bva-p[0]/(1+exp(-p[1]-p[2]*ava));
}
}
double corrcoef(double *d1, double *d2, int dataL)
{
int i;
double xy=0, x=0, y=0, xsum=0, ysum=0;
double corrc;
for (i=0; i<dataL; i++)
{
xsum += d1;
ysum += d2;
}
xsum /= dataL;
ysum /= dataL;
for (i=0; i<dataL; i++)
{
d1 -= xsum;
d2 -= ysum;
x += d1 * d1;
y += d2 * d2;
xy += d1 * d2;
}
corrc = abs(xy) / (sqrt(x) * sqrt(y));
return corrc;
}
int main()
{
register int i;
int ret;
double p[5], // 5 is max(2, 3, 5)
x[16]; // 16 is max(2, 3, 5, 6, 16)
int m, n;
double opts[LM_OPTS_SZ], info[LM_INFO_SZ];
char *probname[]={
"Rosenbrock function",
"modified Rosenbrock problem",
"Powell's function",
"Wood's function",
"Meyer's (reformulated) problem",
"Osborne's problem",
"helical valley function",
"Boggs & Tolle's problem #3",
"Hock - Schittkowski problem #28",
"Hock - Schittkowski problem #48",
"Hock - Schittkowski problem #51",
"Hock - Schittkowski problem #01",
"Hock - Schittkowski modified problem #21",
"hatfldb problem",
"hatfldc problem",
"equilibrium combustion problem",
"Hock - Schittkowski modified #1 problem #52",
"Schittkowski modified problem #235",
"Boggs & Tolle modified problem #7",
"Hock - Schittkowski modified #2 problem #52",
"Hock - Schittkowski modified problem #76",
};
opts[0]=LM_INIT_MU; opts[1]=1E-15; opts[2]=1E-15; opts[3]=1E-20;
opts[4]= LM_DIFF_DELTA;
// char *reason[]={
// "",
// "stopped by small gradient J^T e",
// "stopped by small Dp",
// "stopped by itmax",
// "singular matrix. Restart from current p with increased mu",
// "no further error reduction is possible. Restart with increased mu",
// "stopped by small ||e||_2"
// };
//
// int i, m, n, ret;
// n = 2; m = 2;
// double p[10], x[10];
// double info[100];
//
// char JOBU;
// char JOBVT;
//
//// int i;
//
// //数据类型integer是fortran里的。这里在C++下可以使用的原因是f2c.h文件中已经作了定义
// integer M = SIZE;
// integer N = SIZE;
// integer LDA = M;
// integer LDU = M;
// integer LDVT = N;
// integer LWORK;
// integer INFO;
//
// integer mn = min( M, N );
//
// integer MN = max( M, N );
//
// double a[SIZE*SIZE] = { 16.0, 5.0, 9.0 , 4.0, 2.0, 11.0, 7.0 , 14.0, 3.0, 10.0, 6.0, 15.0, 13.0, 8.0, 12.0, 1.0};
// double s[SIZE];
// double wk[201];
// double uu[SIZE*SIZE];
// double vt[SIZE*SIZE];
//
// JOBU = 'A';
//
// JOBVT = 'A';
//
// LWORK = 201;
/* Subroutine int dgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
doublereal *a, integer *lda, doublereal *s, doublereal *u, integer *
ldu, doublereal *vt, integer *ldvt, doublereal *work, integer *lwork,
integer *info)
*/
//dgesvd_( &JOBU, &JOBVT, &M, &N, a, &LDA, s, uu, &LDU, vt, &LDVT, wk, &LWORK, &INFO);
//printf("INFO=%d \n", INFO );
//for ( i= 0; i< SIZE; i++ ) {
// printf("s[ %d ] = %f\n", i, s[ i ] );
//}
//
//p[0] = 2; p[1] = 3; // initial values
//for(i = 0; i < n; i++) x = 0.0;
//ret=dlevmar_dif(func, p, x, m, n, 1000, NULL, info, NULL, NULL, NULL); // with numerical jacobian
//printf("\nUnconstrained, with numerical jacobian\n--------------------------------------");
//printf("\nLevenberg-Marquardt returned in %d iteration(s): %s\nSolution: ", ret, reason[int(info[6])]);
//for(i=0; i<m; ++i) printf("%.7g ", p);
//printf("\n");
///* Woods's function */
//m=4; n=6;
//p[0]=-3.0; p[1]=-1.0; p[2]=-3.0; p[3]=-1.0;
//for(i=0; i<n; i++) x=0.0;
////x[1]=1.0;
//ret=dlevmar_dif(wood, p, x, m, n, 1000, NULL, info, NULL, NULL, NULL); // no jacobian
//for(i=0; i<m; ++i) printf("%.7g ", p);
//printf("\n");
//m=2; n=2;
//p[0]=2.0; p[1]=1.0;
//for(i=0; i<n; i++) x=0.0;
//ret=dlevmar_dif(ros, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
////ret=dlevmar_dif(ros, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // no Jacobian
//for(i=0; i<m; ++i) printf("%.7g ", p);
//printf("\n");
//double ava[3]={1,2,3};
//double bva[3]={3,6,10};
//for(i=0;i<3;i++)
//b=2/(1+1*exp(-2*a));
//for(i=0; i<m; ++i) printf("%.7g ", b);printf("\n");
m=4; n=4;
p[0]=23;p[1]=23;p[2]=0.0008;p[3]=2;
for(i=0; i<n; i++) x=0.0;
//for(i=0;i<6;i++)
// bva=5/(1+exp(-1-1*ava));
//ret=dlevmar_dif(funa,p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
//ret=dlevmar_dif(Logistic,p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
ret=dlevmar_dif(Weibull,p, x, m, n, 10000, opts, info, NULL, NULL, NULL);
//ret=dlevmar_dif(ros, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // no Jacobian
for(i=0; i<m; ++i) printf("%.7g ", p);
printf("\n");
for(i=0; i<n; ++i) printf("%.7g ", bva);printf("\n");
double cva[6];
for(i=0;i<n;i++)
// cva=p[0]/(1+exp(-p[1]-p[2]*ava));
//cva=p[0]/(1+p[1]*exp(-p[2]*ava)); //Logistic
cva=p[0]-p[1]*exp(-p[2]*pow(ava,p[3]));
for(i=0; i<n; ++i) printf("%.7g ", cva);printf("\n");
double sva=0;
sva = corrcoef(cva,bva,n);
printf("%.7g ", sva);
system("pause");
return 0;
} |
|