numerical recipes（C语言*fp = *dp如何理解）

本文目录

• C语言*fp = *dp如何理解
• 用matlab能求spearman秩相关系数吗
• 有谁用C语言实现过correl函数

C语言*fp = *dp如何理解

Numerical Recipes in C
一书及所附程序，有完整程序。不过我封装了它的C++版本，可以对但参数或多参数求极值，完整的头文件为：
#ifndef __OPTIMIZATION_H__
#define __OPTIMIZATION_H__
//////////////////////////////////////////////////////////////////////////
// class TOptimization
//
// $ 求函数一个或多个参数的最小值
//
// 该类默认对一个参数优化，一般只要输入优化参数数目
// 优化目标函数就可以使用。
//
// ...主要代码来自：
// Numerical Recipes in C++
// The Art of Scientific Computing
// Second Edition
// William H. Press Saul A. Teukolsky
// William T. Vetterling Brian P. Flannery
//
// 中译本：
// C++ 数值算法(第二版) 胡健伟 赵志勇 薛运华 等译
// 电子工业出版社 北京 (2005)
//
// Author: Jian Feng
// Email: fengj@tom.com
// Dec. 9, 2006
//
//////////////////////////////////////////////////////////////////////////
//
// 输入函数：
//
// @MaxIterationStep: 最大迭代次数， 默认 1000
// @ParameterNumbers: 优化参数数目， 默认 1
// @InitMatrix: 初始化矩阵参数(N*N)， 默认
// @Tolerance: 容许公差， 默认 1E-7
//
// 执行函数：
//
// @ExecutePowell: 利用powell方法进行多参数优化
// @ExecuteBrent: 利用brent方法进行单参数优化
//
// 输出函数：
//
// @OptimizatedParameters: 优化结果数据
// @ObjectiveFunctionValue: 目标函数在优化值处的值
//
// 使用示例：
//
// 1. 单参数
// double objfun(double a){
// double sum = 0;
// for(int i = 0; i 《 DataPoints; ++i)
// sum += SQR(Exps - Theo(a));
// }
// double value
// TOptimization opt;
// if(opt.ExecuteBrent(objfun, -10, -1)) opt.OptimizatedParameters(&value);
//
// 2. 多参数
// double objfun(double *a){
// double sum = 0;
// for(int i = 0; i 《 DataPoints; ++i)
// sum += SQR(Exps - Theo(a));
// }
// double value
// TOptimization opt(3);
// double ival = ;
// if(opt.ExecutePowell(objfun, ival)) opt.OptimizatedParameters(value);
//
namespace{
static int ncom; //公用变量
static double *pcom_p; //公用变量
static double *xicom_p; //公用变量
static double (*nrfunc)(double*); //公用函数指针
}
class TOptimization
{
private:
typedef double (*Reff)(double *);
typedef double (*Ptrf)(double );
public:
TOptimization(int n = 1);
~TOptimization()
//主要方法
void ParameterNumbers(int n)
//利用powell方法对一个或多个参数优化
bool ExecutePowell(Reff obj, double *a = 0);
//利用brent方法对一个参数优化，需给出参数所在的区间
bool ExecuteBrent(Ptrf obj, double vFrom = 0, double vTo = 1);
void OptimizatedParameters(double *a)
void OptimizatedParameters(double &a)
//void OptimizatedParameters(double *a){
// if(method) for(int i=0; i《num; ++i) a;
// else *a = vmin;
//}
//其它方法
void InitMatrix(double **m)
{
for(int i=0; i《num; ++i)
for(int j = 0; j《num; ++j)
matx;
setm = true;
}
void MaxIterationStep(int s)
void Tolerance(double eps)
double ObjectiveFunctionValue()const
private:
double brent(double ax, double bx, double cx, Ptrf f, double tol, double &xmin, int &flag);
void mnbrak(double &ax, double &bx, double &cx, double &fa, double &fb, double &fc, Ptrf func);
void linmin(double *p, double *xi, double &fret, Reff func);
bool powell(double *p, double **xi, double ftol, int &iter, double &fret, Reff func);

void shft2(double &a, double &b, const double c)
void shft3(double &a, double &b, double &c, const double d)

double SQR(double x)
void SWAP(double &a, double &b)

double SIGN(const double &a, const double &b)
double MAX(const double &a, const double &b)
void AllocateMemory();
void FreeMemory();
static double f1dim(double x)
{
int j;
double *xt = new double ;
//Vec_Dp &pcom=*pcom_p,&xicom=*xicom_p;
double *pcom = pcom_p, *xicom = xicom_p;
for (j=0;j《ncom;j++)
xt;
//delete xt;
double val = nrfunc(xt);
delete xt;
return val;
}
bool setm; //是否设置优化方向初始矩阵
int num; //优化参数
int ITMAX; //最大迭代数
int iter; //实际迭代步数
int method; //优化方法 0: 1-D brent, 2: M-D Powell
double vmin; //一维优化参数
double ftol; //容许差
double fret; //目标函数值
double *coef; //多维优化参数值
double **matx; //多维优化参数方向的初始值
};
//////////////////////////////////////////////////////////////////////////
inline TOptimization::TOptimization(int n )
{
num = n;
ftol = 1e-7;
ITMAX = 1000;
iter = 0;
fret = 0.;
vmin = 0.;
method = 0;
setm = false;
AllocateMemory();
}
inline void TOptimization::AllocateMemory()
{
pcom_p = new double ;
xicom_p = new double ;
coef = new double ;
matx = new double *;
for(int i = 0; i 《 num; ++i)
{
coef = 0.;
matx;
for(int j = 0; j 《 num; ++j)
matx=(i == j ? 1.0 : 0.0);
}
}
inline void TOptimization::FreeMemory()
{
for(int i = 0; i 《 num; ++i)
{
delete ;
}
delete matx;
delete pcom_p;
delete xicom_p;
delete coef;
}
inline bool TOptimization::ExecutePowell(Reff obj, double *a)
{
method = 1;
if(a)
for(int i = 0; i 《 num; ++i) coef;
return powell(coef, matx, ftol, iter, fret, obj);
}
inline bool TOptimization::ExecuteBrent(Ptrf obj, double vFrom, double vTo)
{
method = 0;
int flag;
double cx, fa, fb, fc;
mnbrak(vFrom,vTo,cx,fa,fb,fc,obj);
fret = brent(vFrom,vTo,cx,obj, ftol,vmin, flag);
return flag ? true : false;
}
inline void TOptimization::mnbrak(double &ax, double &bx, double &cx, double &fa,
double &fb, double &fc, Ptrf func)
{
const double GOLD=1.618034,GLIMIT=100.0,TINY=1.0e-20;
double ulim,u,r,q,fu;
fa=func(ax);
fb=func(bx);
if (fb 》 fa) {
SWAP(ax,bx);
SWAP(fb,fa);
}
cx=bx+GOLD*(bx-ax);
fc=func(cx);
while (fb 》 fc) {
r=(bx-ax)*(fb-fc);
q=(bx-cx)*(fb-fa);
u=bx-((bx-cx)*q-(bx-ax)*r)/
(2.0*SIGN(MAX(fabs(q-r),TINY),q-r));
ulim=bx+GLIMIT*(cx-bx);
if ((bx-u)*(u-cx) 》 0.0) {
fu=func(u);
if (fu 《 fc) {
ax=bx;
bx=u;
fa=fb;
fb=fu;
return;
} else if (fu 》 fb) {
cx=u;
fc=fu;
return;
}
u=cx+GOLD*(cx-bx);
fu=func(u);
} else if ((cx-u)*(u-ulim) 》 0.0) {
fu=func(u);
if (fu 《 fc) {
shft3(bx,cx,u,cx+GOLD*(cx-bx));
shft3(fb,fc,fu,func(u));
}
} else if ((u-ulim)*(ulim-cx) 》= 0.0) {
u=ulim;
fu=func(u);
} else {
u=cx+GOLD*(cx-bx);
fu=func(u);
}
shft3(ax,bx,cx,u);
shft3(fa,fb,fc,fu);
}
}
inline double TOptimization::brent(double ax, double bx, double cx,
Ptrf f, double tol, double &xmin, int &flag)
{
flag = 1;
const double CGOLD=0.3819660;
const double ZEPS=1.0e-20;
int iter;
double a,b,d=0.0,etemp,fu,fv,fw,fx;
double p,q,r,tol1,tol2,u,v,w,x,xm;
double e=0.0;
a=(ax 《 cx ? ax : cx);
b=(ax 》 cx ? ax : cx);
x=w=v=bx;
fw=fv=fx=f(x);
for (iter=0;iter《ITMAX;iter++) {
xm=0.5*(a+b);
tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
if (fabs(x-xm) 《= (tol2-0.5*(b-a))) {
xmin=x;
return fx;
}
if (fabs(e) 》 tol1) {
r=(x-w)*(fx-fv);
q=(x-v)*(fx-fw);
p=(x-v)*q-(x-w)*r;
q=2.0*(q-r);
if (q 》 0.0) p = -p;
q=fabs(q);
etemp=e;
e=d;
if (fabs(p) 》= fabs(0.5*q*etemp) || p 《= q*(a-x) || p 》= q*(b-x))
d=CGOLD*(e=(x 》= xm ? a-x : b-x));
else {
d=p/q;
u=x+d;
if (u-a 《 tol2 || b-u 《 tol2)
d=SIGN(tol1,xm-x);
}
} else {
d=CGOLD*(e=(x 》= xm ? a-x : b-x));
}
u=(fabs(d) 》= tol1 ? x+d : x+SIGN(tol1,d));
fu=f(u);
if (fu 《= fx) {
if (u 》= x) a=x; else b=x;
shft3(v,w,x,u);
shft3(fv,fw,fx,fu);
} else {
if (u 《 x) a=u; else b=u;
if (fu 《= fw || w == x) {
v=w;
w=u;
fv=fw;
fw=fu;
} else if (fu 《= fv || v == x || v == w) {
v=u;
fv=fu;
}
}
}
flag = 0;
xmin=x;
return fx;
}
inline void TOptimization::linmin(double *p, double *xi, double &fret, Reff func)
{
int j, flag;
const double TOL=1.0e-8;
double xx,xmin,fx,fb,fa,bx,ax;
int n=num;
ncom=n;
//pcom_p=new Vec_Dp(n);
//xicom_p=new Vec_Dp(n);
nrfunc=func;
//Vec_Dp &pcom=*pcom_p,&xicom=*xicom_p;
double *pcom = pcom_p, *xicom = xicom_p;
for (j=0;j《n;j++) {
pcom;
xicom;
}
ax=0.0;
xx=1.0;
mnbrak(ax,xx,bx,fa,fx,fb,f1dim);
fret=brent(ax,xx,bx,f1dim,TOL,xmin, flag);
for (j=0;j《n;j++) {
xi *= xmin;
p;
}
//delete xicom_p;
//delete pcom_p;
}
inline bool TOptimization::powell(double *p, double **xi, double ftol, int &iter,
double &fret, Reff func)
{
const int ITMAX=500;
const double TINY=1.0e-20;
int i,j,ibig;
double del,fp,fptt,t;
int n=num;
//Vec_Dp pt(n),ptt(n),xit(n);
double *pt, *ptt, *xit;
for(i = 0; i 《 n; ++i)
{
pt = new double ;
ptt = new double ;
xit = new double ;
}
fret=func(p);
for (j=0;j《n;j++) pt;
for (iter=0;;++iter) {
fp=fret;
ibig=0;
del=0.0;
for (i=0;i《n;i++) {
for (j=0;j《n;j++) xit;
fptt=fret;
linmin(p,xit,fret,func);
if (fptt-fret 》 del) {
del=fptt-fret;
ibig=i+1;
}
}
if (2.0*(fp-fret) 《= ftol*(fabs(fp)+fabs(fret))+TINY) {
delete pt;
delete ptt;
delete xit;
return true;
}
if (iter == ITMAX)
{
delete pt;
delete ptt;
delete xit;
return false;
//cerr《《"powell exceeding maximum iterations.";
}
for (j=0;j《n;j++) {
ptt;
xit;
pt;
}
fptt=func(ptt);
if (fptt 《 fp) {
t=2.0*(fp-2.0*fret+fptt)*SQR(fp-fret-del)-del*SQR(fp-fptt);
if (t 《 0.0) {
linmin(p,xit,fret,func);
for (j=0;j《n;j++) {
xi;
xi;
}
}
}
}
}
#endif

用matlab能求spearman秩相关系数吗

老外的例子，不需要看懂英文，有算例。
Statistics Toolbox
***隐藏网址***
Example:
》》 x = ’;
》》 y = ’;
》》 = spear(x,y)
》》 =spear(x,y)
r =
0.8250 -0.6000
t =
2.5285 -1.2990
p =
0.0855 0.2848

有谁用C语言实现过correl函数

很经典的 Numerical Recipes in C  里有。网上可以查到。

网页链接

例子：

void correl(data1,data2,n,ans)
float data1;
int n;
{
int no2,i;
float dum,*fft,*vector();
void twofft(),realft(),free_vector();

fft=vector(1,2*n);
twofft(data1,data2,fft,ans,n);
no2=n/2;
for (i=2;i《=n+2;i+=2) {
ans)/no2;
ans)/no2;
}
ans;
realft(ans,no2,-1);
free_vector(fft,1,2*n);
}

问题在于还要其他函数代码。你需要有统计数学的基本知识，才能运用和改进满足自己的项目需要。