花如月,谢谢哦,我已经使用C++实现了ksdensity的正态分布密度功能。
此为matlab的ksdensity的C++实现,只考虑了正态分布
// density.h: interface for the Cdensity class.
//
//////////////////////////////////////////////////////////////////////
#if !defined(AFX_DENSITY_H__DF52D4E3_8E9D_45AE_8BD8_97062CCEB182__INCLUDED_)
#define AFX_DENSITY_H__DF52D4E3_8E9D_45AE_8BD8_97062CCEB182__INCLUDED_
#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000
#include "MATRIX.h"
class Cdensity
{
public:
CMatrix normal(CMatrix z);
CMatrix linspace(double ,double ,int m=100);
void ksdensity(CMatrix y,int m,CMatrix &x,CMatrix &f);
Cdensity();
virtual ~Cdensity();
};
#endif // !defined(AFX_DENSITY_H__DF52D4E3_8E9D_45AE_8BD8_97062CCEB182__INCLUDED_)
// density.cpp: implementation of the Cdensity class.
//
//////////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include "density.h"
#include "math.h"
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif
//#define SefeDeleteArray(ps) if(ps!=NULL) { delete []ps;ps=NULL;}
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
Cdensity::Cdensity()
{
}
Cdensity::~Cdensity()
{
}
//y 必须是列向量,x是返回区间,f是返回密度函数
void Cdensity::ksdensity(/*in*/CMatrix y,/*in*/int m,/*out*/CMatrix &x,/*out*/CMatrix &f)
{
//y 必须是列向量
double n;
int size;
if(1==y.GetColumnNum()&&1!=y.GetRowNum())
{
size=y.GetRowNum();
n=y.GetRowNum();
double ymin,ymax;
ymin=y.Min();
ymax=y.Max();
//中值,先排序再取中值
double sig;
double u;
sig=y.median2()/0.6745;
if(sig<=0)
{
sig=ymax-ymin;
}
if(sig>0)
{
u=sig * pow((4/(3*n)),0.2);
}
else
{
u=1;
}
// double m=100; //也可以使用传递数值,随时变更
//产生100个线性数
if(m<=2)
{
m=100;
}
x=linspace(ymin-2*u,ymax+2*u,m);
CMatrix z,z1,z2;
//第一项按行,第二项按列
z1=z1.repmat(size,x,true);
z2=z2.repmat(m,y);
z=(z1-z2)/u;
z1=normal(z);
f=z1.Sum(1);
f=f/(n*u);
}
}
CMatrix Cdensity::linspace(double min,double max,int m)//返回行向量
{
CMatrix result;
result.Init(1,m);
double step;
int i;
step=(max-min)/(m-1);
for (i=0;i<m-1;++i)
{
result.SetElement(0,i,min+i*step);
}
result.SetElement(0,m-1,max);
return result;
}
CMatrix Cdensity::normal(CMatrix z)
{
int row,col;
row=z.GetRowNum();
col=z.GetColumnNum();
int i,j;
double ftemp;
CMatrix f;
f.Init(row,col);
for(i=0;i<row;++i)
{
for (j=0;j<col;++j)
{
ftemp = exp(-0.5 * z.GetElement(i,j)*z.GetElement(i,j))/ sqrt(2*3.1415926);
f.SetElement(i,j,ftemp);
}
}
return f;
}
CMatrix CMatrix::repmat(int m,CMatrix A,bool flag)
{
CMatrix result;
int row,col,i,j;
row=A.GetRowNum();
col=A.GetColumnNum();
result.Init();
if (flag)//表示为行向量
{
if(1!=row)
{
return result;
}
result.Init(m,col);
for (i=0;i<m;++i)
{
for(j=0;j<col;j++)
{
result.SetElement(i,j,A.GetElement(0,j));
}
}
}
else//表示为列向量
{
if(1!=col)
{
return result;
}
result.Init(row,m);
for (i=0;i<m;++i)
{
for(j=0;j<row;j++)
{
result.SetElement(j,i,A.GetElement(j,0));
}
}
}
return result;
}
本人在求小样本数据的概率密度曲线,困扰了好几天,前天在花如月建议下使用ksdensity函数,在仔细分析了ksdensity的代码后,改用C++实现。以上为核心代码。
[ 本帖最后由 ChaChing 于 2010-1-24 11:18 编辑 ] |