之前遇到的一个客户需求就是计算标准化降水指数(SPI)。

主要参考的论文为：洪兴骏，等 “标准化降水指数SPI分布函数的适用性研究”。原文可以很容易搜索到。

论文将计算方法其实写的比较清楚了，本文主要提供其C++实现。

头文件

//SPI_Calculator.h

#include

#include

#include

#include

#include

using namespace std;

#ifndef SPICALCULATOR_H

#define SPICALCULATOR_H

class SPI_Calculator{

public:

/**

path is the file that provides the data of rain fall

file format:

yearamount of precipitation

yearamount of precipitation

.

.

.

*/

SPI_Calculator(char* path);

/**

parameter is the year to get the SPI

return SPI

*/

double getSPI(int year);

double f_probabilityDensity(double x);

double f_gamma(double z, unsigned Ntermsx);

double integral(double a,double b);

void printYears();

void printSPI();

void test();

private:

bool readData(char* path);

void calculateParameter();

vector> dataTable;

double dataCounter,zeroCounter;

double xAver; //average number of samples

double c0,c1,c2,d1,d2,d3; //const parameter

double para_A, para_Gamma, para_Beta;

};

#endif

实现文件：

#include"SPI_Calculator.h"

SPI_Calculator::SPI_Calculator(char* path){

c0 = 2.515517;

c1 = 0.802853;

c2 = 0.010328;

d1 = 1.432788;

d2 = 0.189269;

d3 = 0.001308;

dataCounter = 0;

zeroCounter = 0;

if(!readData(path)){

cout<

system("pause");

exit(0);

}

calculateParameter();

}

bool SPI_Calculator::readData(char* path){

fstream dataIn(path,ios::in);

if(!dataIn.is_open()){

cout<

return false;

}

pair tempPair;

while(!dataIn.eof()){

dataIn>>tempPair.first;

dataIn>>tempPair.second;

dataTable.push_back(tempPair);

dataCounter ++;

if(tempPair.second==0)

zeroCounter ++;

if(tempPair.second < 0){

cout<

dataIn.close();

return false;

}

}

dataIn.close();

return true;

}

void SPI_Calculator::printYears(){

cout<

for(int i = 0; i < dataCounter; i ++)

cout<

cout<

}

void SPI_Calculator::printSPI(){

cout<

for(int i = 0; i < dataCounter; i ++){

cout<

}

cout<

}

void SPI_Calculator::calculateParameter(){

double sum = 0;

for(int i = 0; i < dataCounter; i ++)

sum += dataTable[i].second;

xAver = sum / (dataCounter-zeroCounter);

para_A = log(xAver);

double temp = 0.0;

for(int i = 0; i < dataCounter; i ++)

if((dataTable[i].second!=0))

temp += log(dataTable[i].second);

para_A -= (temp / (dataCounter-zeroCounter));

para_Gamma = (1 + sqrt(1+4*para_A/3))/ (4* para_A);

para_Beta = xAver / para_Gamma;

}

double SPI_Calculator::getSPI(int year){

double perciption;

int i = 0;

for(; i < dataCounter; i ++)

if(year == dataTable[i].first){

perciption = dataTable[i].second;

break;

}

if(i == dataCounter){

cout<

return -999;

}

double F = 0.0;

if(perciption == 0){

F=(zeroCounter / dataCounter);

}

else

F = integral(0.0000001, perciption);

double S = 0.0;

if(F > 0.5){

F = 1-F;

S = 1;

}

else

S = -1;

double t = sqrt(log(1/(F*F)));

double Z = S*(t - (c0+c1*t+c2*t*t)/(1+d1*t+d2*t*t+d3*t*t*t));

return Z;

}

double SPI_Calculator::f_probabilityDensity(double x){

double E = 2.71828;

double t_a = 1/(pow(para_Beta,para_Gamma)* f_gamma(para_Gamma,1000));

double t_b = pow(x, para_Gamma-1);

double t_c = pow(E, -x/para_Beta);

return (t_a *t_b *t_c);

}

double SPI_Calculator::f_gamma(double z, unsigned Nterms){

double g = 0.577216;// Euler-Mascheroni constant

double retVal = z*exp(g*z);

// apply products

for(unsigned n=1; n<=Nterms; ++n)

retVal *= (1 + z/n)*exp(-z/n);

retVal = 1.0/retVal;// invert

return retVal;

}

double SPI_Calculator::integral(double a,double b) {

int N = 100000*para_A*para_A;

double s = 0,h;

s=(f_probabilityDensity(a) + f_probabilityDensity(b))/2.0;

h=(b-a)/N;

for(int i=1; i

s += f_probabilityDensity(a+i*h);

return (s*h);

}

其实数学的东西实现起来没有特别复杂，两个难点在于积分(概率密度)以及 gamma 函数，在代码里也写清楚了。

使用的时候只需要按代码注释中说明的格式传入每年的降水量就行了。调用 getSPI (int year) 函数可以获取当年的 SPI。

注意，本程序基于合理的数据进行检测。不合理的数据(如每年的降水量完全相同，被视为不合理的数据)将引起意料之外的结果。