Levenberg-Marquardt算法问题 ,请教关于levmar-2.5中dlevmar_dif()的使用
这个是Levenberg-Marquardt算法,找了很多资料,都没有发现他要求提供的函数关系的原型是如何得出的int dlevmar_dif(
void (*func)(double *p, double *hx, int m, int n, void *adata),
double *p, double *x, int m, int n, int itmax, double *opts,
double *info, double *work, double *covar, void *adata);
/* Secant version of the LEVMAR_DER() function above: the Jacobian is approximated with
* the aid of finite differences (forward or central, see the comment for the opts argument)
*/
int LEVMAR_DIF(
void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata), /* functional relation describing measurements. A p \in R^m yields a \hat{x} \inR^n */
LM_REAL *p, /* I/O: initial parameter estimates. On output has the estimated solution */
LM_REAL *x, /* I: measurement vector. NULL implies a zero vector */
int m, /* I: parameter vector dimension (i.e. #unknowns) */
int n, /* I: measurement vector dimension */
int itmax, /* I: maximum number of iterations */
LM_REAL opts, /* I: opts = minim. options [\mu, \epsilon1, \epsilon2, \epsilon3, \delta]. Respectively the
* scale factor for initial \mu, stopping thresholds for ||J^T e||_inf, ||Dp||_2 and ||e||_2 and
* the step used in difference approximation to the Jacobian. Set to NULL for defaults to be used.
* If \delta<0, the Jacobian is approximated with central differences which are more accurate
* (but slower!) compared to the forward differences employed by default.
*/
LM_REAL info,
/* O: information regarding the minimization. Set to NULL if don't care
* info= ||e||_2 at initial p.
* info=[ ||e||_2, ||J^T e||_inf,||Dp||_2, mu/max_ii ], all computed at estimated p.
* info= # iterations,
* info=reason for terminating: 1 - stopped by small gradient J^T e
* 2 - stopped by small Dp
* 3 - stopped by itmax
* 4 - singular matrix. Restart from current p with increased mu
* 5 - no further error reduction is possible. Restart with increased mu
* 6 - stopped by small ||e||_2
* 7 - stopped by invalid (i.e. NaN or Inf) "func" values. This is a user error
* info= # function evaluations
* info= # Jacobian evaluations
* info= # linear systems solved, i.e. # attempts for reducing error
*/
LM_REAL *work, /* working memory at least LM_DIF_WORKSZ() reals large, allocated if NULL */
LM_REAL *covar, /* O: Covariance matrix corresponding to LS solution; mxm. Set to NULL if not needed. */
void *adata) /* pointer to possibly additional data, passed uninterpreted to func.
* Set to NULL if not needed
*/
以上是函数说明,需要的是用这个函数拟合曲线y=c/(1+exp(a+bx))+d;(提供四个点)
以下是在Nemerical Recipes in C这本书中函数所需要的原型,不知道如何转化为dlevmar_dif所要求的形式
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*x;
ex=exp(fac);
arg=(1+ex)*(1+ex);
*y=a/(1+ex)+a;
dyda=(-a)*ex/arg;
dyda=(-a)*x*ex/arg;
dyda=1/(1+ex);
dyda=1;
}
}
那位大侠知道如何使用或者有更好的曲线拟合库的请指教
它的网站 http://www.ics.forth.gr/~lourakis/levmar/
附上安装方法
levmar-2.5安装步骤，测试通过
转载自http://blog.csdn.net/lkbwx/archive/2010/02/19/5311802.aspx
点击http://www.netlib.org/clapack/CLAPACK-3.1.1-VisualStudio.zip 下载clapack
点击http://www.cmake.org/files/v2.8/cmake-2.8.0-win32-x86.exe 下载cmake
点击http://www.ics.forth.gr/~lourakis/levmar/levmar-2.5.tgz 下载LEVMAR
打开CLAPACK-3.1.1-VisualStudio文件夹，先把\LIB\Win32的lib都删了,以免混淆
双击clapack.vcproj打开工程项目文件，下面的各编译步骤都编译成debug模式
依次编译F2CLIBS,tmglib,blas,CLAPACK,结果是在\LIB\Win32下依次生成了libf2cd.lib,tmglibd.lib,BLASd.lib,clapackd.lib
使用cmake作用于levmar-2.5文件夹，得到工程项目文件LEVMAR.vcproj，打开它
工具->选项->项目和解决方案->vc++目录->包含文件处添加 \INCLUDE目录
工具->选项->项目和解决方案->vc++目录->库文件处添加 \LIB\Win32目录
项目属性->连接器->输入->附加依赖项libf2cd.lib BLASd.lib clapackd.lib tmglibd.lib
项目属性->连接器->输入->忽略特定库Libcmtd.lib
欲在LEVMAR中调用CLAPACK，还需要打开LEVMAR的misc_core.c文件，然后在文件起始处写上#include "blaswrap.h"
参考
http://www.cnblogs.com/random_walk/archive/2009/11/13/1602682.html
http://www.wilmott.com/messageview.cfm?catid=34&threadid=51193 找到了例子中的原型，发现其实很简单，但是初值的选择还是有一定的难度，网上资料很少，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={1,2,3,4,5,6};
double bva={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)*(1.0-p) + ROSD*(p-p*p)*(p-p*p));
//}
//
//void jacros(double *p, double *jac, int m, int n, void *data)
//{
// register int i, j;
//
// for(i=j=0; i<n; ++i){
// jac=(-2 + 2*p-4*ROSD*(p-p*p)*p);
// jac=(2*ROSD*(p-p*p));
// }
//}
//
//#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 - 3.1415926) * (p - 3.1415926) +
// (p - 2.0)* (p - 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 - p*p);
// x=1.0 - p;
// x=sqrt(90.0)*(p - p*p);
// x=1.0 - p;
// x=sqrt(10.0)*(p+p - 2.0);
// x=(p - p)/sqrt(10.0);
// }
//}
//void func(double *p, double *x, int m, int n, void *data)
//{
// int i;
// double a={1,2,3};
// double b;
// 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*a*a*a+p*a*a+p*a);
// x =b-(p*a*a*a+p*a*a+p*a);
// x =b-(p*a*a*a+p*a*a+p*a);
// }
//}
//void jacfunc(double *p, double *jac, int m, int n, void *data)
//{
// register int i, j;
//
// for(i=j=0; i<n; ++i){
// jac=-1;
// jac=-1;
//
// jac=-4;
// jac=-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*x;
// ex=exp(fac);
// arg=(1+ex)*(1+ex);
// *y=a/(1+ex)+a;
// dyda=(-a)*ex/arg;
// dyda=(-a)*x*ex/arg;
// dyda=1/(1+ex);
// dyda=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);
x = bva-(p-p*exp(-p*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=-1;
//// jac=exp(-p*pow(ava,p));
//// jac=-p*exp(-p*pow(ava,p))*-p*pow(ava,p)*(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/(1+p*exp(-p*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/(1+exp(-p-p*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;
intret;
double p, // 5 is max(2, 3, 5)
x; // 16 is max(2, 3, 5, 6, 16)
int m, n;
double opts, info;
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=LM_INIT_MU; opts=1E-15; opts=1E-15; opts=1E-20;
opts= 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, x;
// double info;
//
// 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 = { 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;
// double wk;
// double uu;
// double vt;
//
// 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 = 2; p = 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)]);
//for(i=0; i<m; ++i) printf("%.7g ", p);
//printf("\n");
///* Woods's function */
//m=4; n=6;
//p=-3.0; p=-1.0; p=-3.0; p=-1.0;
//for(i=0; i<n; i++) x=0.0;
////x=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=2.0; p=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={1,2,3};
//double bva={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=23;p=23;p=0.0008;p=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;
for(i=0;i<n;i++)
// cva=p/(1+exp(-p-p*ava));
//cva=p/(1+p*exp(-p*ava));//Logistic
cva=p-p*exp(-p*pow(ava,p));
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;
}
页:
[1]