彭汉川版mRMR特征选择算法网址:https://www.mathworks.com/matlabcentral/fileexchange/14608-mrmr-feature-selection-using-mutual-information-computation
estpa函数是为了计算单个向量a的概率分布,其大概思想是:
当a=[1 2 1 2 1]时,其概率分布为 p(x=1)=3/5=0.6, p(x=2)=2/5=0.4,所得P=[0.6 0.4]
当a=[1 5 1 2 1]时,其概率分布为p(x=1)=3/5=0.6, p(x=2)=1/5=0.2, p(x=3)=0/5=0, p(x=4)=0/5=0, p(x=5)=1/5=0.2, 所得P=[0.6 0.2 0 0 0.2]
//=========================================================
//
//This is a prog in the MutualInfo 0.9 package written by
// Hanchuan Peng.
//
//Disclaimer: The author of program is Hanchuan Peng
// at <penghanchuan@yahoo.com> and <phc@cbmv.jhu.edu>.
//
//The CopyRight is reserved by the author.
//
//Last modification: April/19/2002
//
//========================================================
//
// estpa.cpp
// Calculate the histgram/probability of one vector/image
// By Hanchuan Peng
// Dec/2000
// April/8/2002
#include "miinclude.h"
//return the number of states
template <class T> void copyvecdata(T * srcdata, long len, int * desdata, int& nstate, int &minn, int& maxx);
template <class T> void copyvecdata(T * srcdata, long len, int * desdata, int& nstate, int &minn, int& maxx)
{
//输入输出变量定义:
//srcdata:输入向量
//len:输入向量的长度
//desdata:输出向量地址
//nstate:最大值-最小值+1
//例如输入矩阵为a=[1 2 1 2 1]
if(!srcdata || !desdata)
{
printf("NULL points in copyvecdata()!\n");
return;
}
long i;
//note: originally I added 0.5 before rounding, however seems the negative numbers and
// positive numbers are all rounded towarded 0; hence int(-1+0.5)=0 and int(1+0.5)=1;
// This is unwanted because I need the above to be -1 and 1.
// for this reason I just round with 0.5 adjustment for positive and negative differently
//copy data
//int minn,maxx;
if (srcdata[0]>0)
maxx = minn = int(srcdata[0]+0.5);
else
maxx = minn = int(srcdata[0]-0.5);
int tmp;
double tmp1;
for (i=0;i<len;i++)
{
tmp1 = double(srcdata[i]);//i=0时,tmp1=1,tmp=2;i=1时,tmp1=2,tmp=3;以此类推
tmp = (tmp1>0)?(int)(tmp1+0.5):(int)(tmp1-0.5);//round to integers将原始值+/-0.5再取整
minn = (minn<tmp)?minn:tmp;
maxx = (maxx>tmp)?maxx:tmp;
desdata[i] = tmp;
}//得到minn=2,maxx=3,desdata[]=[2 3 2 3 2]
//make the vector data begin from 0 (i.e. 1st state)
for (i=0;i<len;i++)
{
desdata[i] -= minn;
}//desdata[]=[0 1 0 1 0]
//return the #state
nstate = (maxx-minn+1);//nstate=3-2+1=2
return;
}
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])//matlab调用主函数
{
// check the arguments检查输入参数是否正确
//其中nlhs为输出参数个数
//plhs[]为输出参数数组
//nrhs为输入参数个数
//prhs[]为输入参数数组
if(nrhs != 1 && nrhs !=2 && nrhs!=3)//输入参数不是1/2/3时报错
mexErrMsgTxt("Usage [marginprob,statelist,cumsumlist] = progname(vector1, maxstatenum, b_returnprob). \n(Both vectors can be images). Max range handled: INT type of the OS");
if(nlhs > 3)//输出参数>3时报错
mexErrMsgTxt("Too many output argument <marginalprob_list, statelist, cumsumlist>.");
if (!mxIsInt8(prhs[0]) && !mxIsUint8(prhs[0]) && !mxIsDouble(prhs[0]) ) //输入参数数据类型不对
mexErrMsgTxt("The first input argument must be types of INT8 or UINT8 or DOUBLE.");
//get and check size information
long i,j;
//输入数据处理
void *img1 = (void *)mxGetData(prhs[0]);//获取非数字矩阵的数据值
long len1 = mxGetNumberOfElements(prhs[0]);//获取数据的数量
mxClassID type1 = mxGetClassID(prhs[0]);//获取数据的类型
if (!img1 || !len1)
mexErrMsgTxt("The input vector is invalid.");
int b_findstatenum = 1;
int nstate1 = 0;
if (nrhs>=2)
{
b_findstatenum = 0;
long MaxGrayLevel = (long) mxGetScalar(prhs[1]);
nstate1 = MaxGrayLevel;
if (MaxGrayLevel<=1)
{
printf("The argument #state is invalid. This program will decide #state itself.\n");
b_findstatenum = 1;
}
}
int b_returnprob = 1;
if (nrhs>=3)
{
b_returnprob = (mxGetScalar(prhs[2])!=0);
}
//copy data into new INT type array (hence quantization) and then reange them begin from 0 (i.e. state1)
int * vec1 = new int[len1];
int nrealstate1=0, minn,maxx;
switch(type1)
{
case mxINT8_CLASS: copyvecdata((char *)img1,len1,vec1,nrealstate1,minn,maxx); break;
//调用前面的数字取整函数
case mxUINT8_CLASS: copyvecdata((unsigned char *)img1,len1,vec1,nrealstate1,minn,maxx); break;
case mxDOUBLE_CLASS: copyvecdata((double *)img1,len1,vec1,nrealstate1,minn,maxx); break;
}//沿用前面的例子vec1=[0 1 0 1 0],nrealstate1=2
//update the #state when necessary
if (nstate1<nrealstate1)
{
nstate1 = nrealstate1;
//printf("First vector #state = %i\n",nrealstate1);
}
//generate the marginal-distribution list
//生成编辑分布
plhs[0] = mxCreateDoubleMatrix(nstate1,1,mxREAL);
double *ha = (double *) mxGetPr(plhs[0]);
for (i=0; i<nstate1;i++)
{
ha[i] = 0;
}//所有值清零
for (i=0;i<len1;i++)
{
ha[vec1[i]] += 1;
//ha[vec1[0]]=ha[0]=1;ha[vec1[1]]=ha[1]=1;循环之后ha[0]=3,ha[1]=2
}
//return the probabilities, otherwise return count numbers
if(b_returnprob)
{
for (i=0; i<nstate1;i++)
{
ha[i] /= len1;
//ha[0]=3/5=0.6;ha[1]=2/5=0.4
}
}
//return more information
if (nlhs>=2) //the second return value is the state list
{
plhs[1] = mxCreateDoubleMatrix(nstate1,1,mxREAL);
double * pstate = (double *)mxGetPr(plhs[1]);
for (i=0;i<nstate1;i++)
pstate[i] = minn + i;
}
if (nlhs>=3) //the third return value is the cumsum prob list
{
plhs[2] = mxCreateDoubleMatrix(nstate1,1,mxREAL);
double * pcumsum = (double *)mxGetPr(plhs[2]);
pcumsum[0] = ha[0];
for (i=1;i<nstate1;i++)
pcumsum[i] = pcumsum[i-1] + ha[i];
}
//free memory
if (vec1) delete []vec1;
return;
}