彭汉川版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]
当b=[2 1 2 1 1]时,其概率分布为p(x=1)=3/5=0.6, p(x=2)=1/5=0.4, 所得P=[0.6, 0.4]
则a,b的联合概率分布为p(a=1,b=1)=1/5=0.2, p(a=1,b=2)=2/5=0.4, p(a=2, b=1)=2/5=0.4, p(a=2,b=2)=0,所得P=[0.2, 0.4; 0.4, 0]
//=========================================================
//
//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
//
//========================================================
//
// estpab.cpp
// Calculate the joint histgram/probability of two vectors/images
//计算两个向量的联合概率
// 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);
template <class T> void copyvecdata(T * srcdata, long len, int * desdata, int& nstate)//数据处理函数
{
//输入变量
//srcdata:输入数据地址
//len:输入数据长度
//desdata:输出数据地址--数据数据值减去最小值
//nstate:最大值减去最小值+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]);
tmp = (tmp1>0)?(int)(tmp1+0.5):(int)(tmp1-0.5);//round to integers将值转换为整数值
minn = (minn<tmp)?minn:tmp;//查找最小值
maxx = (maxx>tmp)?maxx:tmp;//查找最大值
desdata[i] = tmp;
// printf("%i ",desdata[i]);
}
//printf("\n");
//make the vector data begin from 0 (i.e. 1st state)
for (i=0;i<len;i++)
{
desdata[i] -= minn;//desdata里面存储的值为实际值减掉最小值
}
//return the #state
nstate = (maxx-minn+1);//nstate里存储的值为最大值减掉最小值
return;
}
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
// check the arguments检查参数
if(nrhs != 3 && nrhs !=2 && nrhs!=4)
mexErrMsgTxt("Usage [jointprob_table, marginprob_1, marginprob2] = progname(vector1, vector2, maxstatenum, b_returnprob). \n(Both vectors can be images). Max range handled: INT type of the OS");
if(nlhs > 3)
mexErrMsgTxt("Too many output argument <jointprob_table>.");
if (!mxIsInt8(prhs[0]) && !mxIsUint8(prhs[0]) && !mxIsDouble(prhs[0]) )
mexErrMsgTxt("The first input argument must be types of INT8 or UINT8 or DOUBLE.");
if (!mxIsInt8(prhs[1]) && !mxIsUint8(prhs[1]) && !mxIsDouble(prhs[1]) )
mexErrMsgTxt("The second 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]);
void *img2 = (void *)mxGetData(prhs[1]);
long len2 = mxGetNumberOfElements(prhs[1]);
mxClassID type2 = mxGetClassID(prhs[1]);
if (!img1 || !img2 || !len1 || !len2)
mexErrMsgTxt("At least one of the input vectors is invalid.");
if (len1!=len2)
mexErrMsgTxt("The two vectors/images should have the same length.");
int b_findstatenum = 1;
int nstate1 = 0, nstate2 = 0;
if (nrhs>=3)
{
b_findstatenum = 0;
long MaxGrayLevel = (long) mxGetScalar(prhs[2]);
nstate1 = nstate2 = 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>=4)
{
b_returnprob = (mxGetScalar(prhs[3])!=0);
}
//copy data into new INT type array (hence quantization) and then reange them begin from 0 (i.e. state1)
//新建两个变量vec1和vec2分别用来处理输入的两个变量
int * vec1 = new int[len1];
int * vec2 = new int[len2];
int nrealstate1=0, nrealstate2=0;
switch(type1)
{
case mxINT8_CLASS: copyvecdata((char *)img1,len1,vec1,nrealstate1); break;//vec1为整数化并且减去最小值的数据
case mxUINT8_CLASS: copyvecdata((unsigned char *)img1,len1,vec1,nrealstate1); break;
case mxDOUBLE_CLASS: copyvecdata((double *)img1,len1,vec1,nrealstate1); break;
}
switch(type2)
{
case mxINT8_CLASS: copyvecdata((char *)img2,len2,vec2,nrealstate2); break;//vec2为整数化并且减去最小值的数据
case mxUINT8_CLASS: copyvecdata((unsigned char *)img2,len2,vec2,nrealstate2); break;
case mxDOUBLE_CLASS: copyvecdata((double *)img2,len2,vec2,nrealstate2); break;
}
//update the #state when necessary
if (nstate1<nrealstate1)
{
nstate1 = nrealstate1;
// printf("First vector #state = %i\n",nrealstate1);
}
if (nstate2<nrealstate2)
{
nstate2 = nrealstate2;
// printf("Second vector #state = %i\n",nrealstate2);
}
//generate the joint-distribution table
//生成联合分布表
plhs[0] = mxCreateDoubleMatrix(nstate1,nstate2,mxREAL);//生成nstate1行,nstate2列的mxREAL矩阵
double *hab = (double *) mxGetPr(plhs[0]);
double **hab2d = new double * [nstate2];
for(j=0;j<nstate2;j++)
hab2d[j] = hab + (long)j*nstate1;
for (i=0; i<nstate1;i++)
for (j=0; j<nstate2;j++)
{
hab2d[j][i] = 0;
}
for (i=0;i<len1;i++)
{
//old method -- slow
// indx = (long)(vec2[i]) * nstate1 + vec1[i];
// hab[indx] += 1;
//new method -- fast
hab2d[vec2[i]][vec1[i]] += 1;
}
//return the probabilities, otherwise return count numbers
if(b_returnprob)
{
for (i=0; i<nstate1;i++)
for (j=0; j<nstate2;j++)
{
hab2d[j][i] /= len1;
}
}
//for nlhs>=2
if (nlhs>=2)
{
plhs[1] = mxCreateDoubleMatrix(nstate1,1,mxREAL);
double *ha = (double *)mxGetPr(plhs[1]);
for (i=0;i<nstate1;i++) {ha[i] = 0;}
for (i=0;i<nstate1;i++)
for (j=0;j<nstate2;j++)
{
ha[i] += hab2d[j][i];
}
}
if (nlhs>=3)
{
plhs[2] = mxCreateDoubleMatrix(nstate2,1,mxREAL);
double *hb = (double *)mxGetPr(plhs[2]);
for (j=0;j<nstate2;j++) {hb[j] = 0;}
for (i=0;i<nstate1;i++)
for (j=0;j<nstate2;j++)
{
hb[j] += hab2d[j][i];
}
}
//finish
if (hab2d) {delete []hab2d;}
if (vec1) delete []vec1;
if (vec2) delete []vec2;
return;
}