最近有个期权项目，计算理论价时需要使用标准正态分布变量的累积概率分布函数，excel中可以通过normsdist函数得到该结果，但是项目不考虑调用excel公式，于是只能用java来实现这个公式。

先是考虑把正态分布的那张表搞到程序中，通过查表的方式，小数点三位后面多出来的值使用公式来计算，代码如下

private static double[][] normdist = {

{0.5,0.504,0.508,0.512,0.516,0.5199,0.5239,0.5279,0.5319,0.5359},

{0.5398,0.5438,0.5478,0.5517,0.5557,0.5596,0.5636,0.5675,0.5714,0.5753},

{0.5793,0.5832,0.5871,0.591,0.5948,0.5987,0.6026,0.6064,0.6103,0.6141},

{0.6179,0.6217,0.6255,0.6293,0.6331,0.6368,0.6406,0.6443,0.648,0.6517},

{0.6554,0.6591,0.6628,0.6664,0.67,0.6736,0.6772,0.6808,0.6844,0.6879},

{0.6915,0.695,0.6985,0.7019,0.7054,0.7088,0.7123,0.7157,0.719,0.7224},

{0.7257,0.7291,0.7324,0.7357,0.7389,0.7422,0.7454,0.7486,0.7517,0.7549},

{0.758,0.7611,0.7642,0.7673,0.7703,0.7734,0.7764,0.7794,0.7823,0.7852},

{0.7881,0.791,0.7939,0.7967,0.7995,0.8023,0.8051,0.8078,0.8106,0.8133},

{0.8159,0.8186,0.8212,0.8238,0.8264,0.8289,0.8315,0.834,0.8365,0.8389},

{0.8413,0.8438,0.8461,0.8485,0.8508,0.8531,0.8554,0.8577,0.8599,0.8621},

{0.8643,0.8665,0.8686,0.8708,0.8729,0.8749,0.877,0.879,0.881,0.883},

{0.8849,0.8869,0.8888,0.8907,0.8925,0.8944,0.8962,0.898,0.8997,0.9015},

{0.9032,0.9049,0.9066,0.9082,0.9099,0.9115,0.9131,0.9147,0.9162,0.9177},

{0.9192,0.9207,0.9222,0.9236,0.9251,0.9265,0.9278,0.9292,0.9306,0.9319},

{0.9332,0.9345,0.9357,0.937,0.9382,0.9394,0.9406,0.9418,0.943,0.9441},

{0.9452,0.9463,0.9474,0.9484,0.9495,0.9505,0.9515,0.9525,0.9535,0.9545},

{0.9554,0.9564,0.9573,0.9582,0.9591,0.9599,0.9608,0.9616,0.9625,0.9633},

{0.9641,0.9648,0.9656,0.9664,0.9671,0.9678,0.9686,0.9693,0.97,0.9706},

{0.9713,0.9719,0.9726,0.9732,0.9738,0.9744,0.975,0.9756,0.9762,0.9767},

{0.9772,0.9778,0.9783,0.9788,0.9793,0.9798,0.9803,0.9808,0.9812,0.9817},

{0.9821,0.9826,0.983,0.9834,0.9838,0.9842,0.9846,0.985,0.9854,0.9857},

{0.9861,0.9864,0.9868,0.9871,0.9874,0.9878,0.9881,0.9884,0.9887,0.989},

{0.9893,0.9896,0.9898,0.9901,0.9904,0.9906,0.9909,0.9911,0.9913,0.9916},

{0.9918,0.992,0.9922,0.9925,0.9927,0.9929,0.9931,0.9932,0.9934,0.9936},

{0.9938,0.994,0.9941,0.9943,0.9945,0.9946,0.9948,0.9949,0.9951,0.9952},

{0.9953,0.9955,0.9956,0.9957,0.9959,0.996,0.9961,0.9962,0.9963,0.9964},

{0.9965,0.9966,0.9967,0.9968,0.9969,0.997,0.9971,0.9972,0.9973,0.9974},

{0.9974,0.9975,0.9976,0.9977,0.9977,0.9978,0.9979,0.9979,0.998,0.9981},

{0.9981,0.9982,0.9982,0.9983,0.9984,0.9984,0.9985,0.9985,0.9986,0.9986},

{0.9987,0.999,0.9993,0.9995,0.9997,0.9998,0.9998,0.9999,0.9999,1},

{0.999032,0.999065,0.999096,0.999126,0.999155,0.999184,0.999211,0.999238,0.999264,0.999289},

{0.999313,0.999336,0.999359,0.999381,0.999402,0.999423,0.999443,0.999462,0.999481,0.999499},

{0.999517,0.999534,0.999550,0.999566,0.999581,0.999596,0.999610,0.999624,0.999638,0.999660},

{0.999663,0.999675,0.999687,0.999698,0.999709,0.999720,0.999730,0.999740,0.999749,0.999760},

{0.999767,0.999776,0.999784,0.999792,0.999800,0.999807,0.999815,0.999822,0.999828,0.999885},

{0.999841,0.999847,0.999853,0.999858,0.999864,0.999869,0.999874,0.999879,0.999883,0.999880},

{0.999892,0.999896,0.999900,0.999904,0.999908,0.999912,0.999915,0.999918,0.999922,0.999926},

{0.999928,0.999931,0.999933,0.999936,0.999938,0.999941,0.999943,0.999946,0.999948,0.999950},

{0.999952,0.999954,0.999956,0.999958,0.999959,0.999961,0.999963,0.999964,0.999966,0.999967},

{0.999968,0.999970,0.999971,0.999972,0.999973,0.999974,0.999975,0.999976,0.999977,0.999978},

{0.999979,0.999980,0.999981,0.999982,0.999983,0.999983,0.999984,0.999985,0.999985,0.999986},

{0.999987,0.999987,0.999988,0.999988,0.999989,0.999989,0.999990,0.999990,0.999991,0.999991},

{0.999991,0.999992,0.999992,0.999930,0.999993,0.999993,0.999993,0.999994,0.999994,0.999994},

{0.999995,0.999995,0.999995,0.999995,0.999996,0.999996,0.999996,1.000000,0.999996,0.999996},

{0.999997,0.999997,0.999997,0.999997,0.999997,0.999997,0.999997,0.999998,0.999998,0.999998},

{0.999998,0.999998,0.999998,0.999998,0.999998,0.999998,0.999998,0.999998,0.999999,0.999999},

{0.999999,0.999999,0.999999,0.999999,0.999999,0.999999,0.999999,0.999999,0.999999,0.999999},

{0.999999,0.999999,0.999999,0.999999,0.999999,0.999999,0.999999,0.999999,0.999999,0.999999},

{1.000000,1.000000,1.000000,1.000000,1.000000,1.000000,1.000000,1.000000,1.000000,1.000000}

};

static DecimalFormat format = new DecimalFormat("#.00");

static{

format.setRoundingMode(RoundingMode.DOWN);

}

public static double NORMSDIST(double x)

{

if(x<0 || x>4.99)

{

return 0;

}

double rx = x;

x = Double.valueOf(format.format(x));

int row = (int)(x*100)%10;

int col = (int)(x*10);

double rtn = normdist[col][row];

double step = 0.00001;

for(double i = x+step ; i <= rx ; i+=step)

{

rtn+=N_(i)*step;

}

return rtn;

}

private static double N_(double x)

{

double rsp = (1/Math.sqrt(2*Math.PI)) * Math.exp((-1)*Math.pow(x, 2)/2);

return rsp;

}

但是以上方法只能保证精确到小数点后四位，对于计算理论价还是不够；后面去网上搜索，有个牛逼的哥们用了几行代码就搞定了，而且精确度达到小数点后7位。

public static double NORMSDIST(double a)

{

double p = 0.2316419;

double b1 = 0.31938153;

double b2 = -0.356563782;

double b3 = 1.781477937;

double b4 = -1.821255978;

double b5 = 1.330274429;

double x = Math.abs(a);

double t = 1/(1+p*x);

double val = 1 - (1/(Math.sqrt(2*Math.PI))  * Math.exp(-1*Math.pow(a, 2)/2)) * (b1*t + b2 * Math.pow(t,2) + b3*Math.pow(t,3) + b4 * Math.pow(t,4) + b5 * Math.pow(t,5) );

if ( a 

{

val = 1- val;

}

return val;

}

累积分布函数cdf (Cumulative

Distribution Function)

背景知识：http://www.lifelaf.com/blog/?p=746

语法

y =

cdf('name',x,A,B)

y =

cdf('name',x,A,B,C)

y =

cdf(pd,x)

y =

cdf(___,'upper')

描述

y =

cdf('name',x,A) 计算某种分布(由'name'定义，如'Normal'正态,

'Poisson'泊松,

'T' t分布…)下，x值处的累计分布，A,B,C等为'name'函数的参数

y =

cdf(pd,x) 直接计算概率分布函数pd(probability

distribution) ，在x处的累计分布，实际上，这里的pd

已被'name',

A定义好，举栗如下：

%

定义一个正态分布函数pd, 均值mu

= 0, 标准差sigma = 1.

mu = 0;

sigma = 1;

pd =

makedist('Normal',mu,sigma);

%

定义x值

x = [-2,-1,0,1,2];

%

计算x值处的累计分布

y = cdf(pd,x)

y

=

0.0228 0.1587 0.5000 0.8413 0.9772

用第一种语句表达相同内容为：

y2 =

0.0228 0.1587 0.5000 0.8413 0.9772

t分布累积分布函数tcdf

(Student'stcumulative

distribution function)

%

事实上就是y = cdf('T',x,A)函数

语法

p =

tcdf(x,nu)p =

tcdf(x,nu,'upper')

描述

计算t分布在x值处的累积分布，nu是t分布的自由度

再举个栗子

mu = 1; % Population mean

sigma = 2;  % Population standard deviation

n = 100; % Sample size

x = normrnd(mu,sigma,n,1);  % Random sample from population

xbar = mean(x);  % Sample mean

s = std(x); % Sample standard deviation

t = (xbar - mu)/(s/sqrt(n)) %这里t分布出现了，正态分布总体与样本均值的差符合t分布

t =

1.0589

p = 1-tcdf(t,n-1) % Probability of larger t-statistic

p =

0.1461

该p值(即t函数的累积分布就是t检验在相同x值处的概率ptest)

[h,ptest] = ttest(x,mu,0.05,'right')

h =

0

ptest =

0.1461

概率密度函数pdf (Probability

density functions)

搞懂了累积分布函数cdf，这个就没什么需要多说了

语法

y =

pdf('name',x,A)

y =

pdf('name',x,A,B)

y =

pdf('name',x,A,B,C)

y =

pdf(pd,x)

举例

%

定义一个正态分布函数pd, 均值mu

= 0, 标准差sigma = 1.

mu =

0;

sigma =

1;

pd =

makedist('Normal',mu,sigma);

%

定义x值

x = [-2

-1 0 1 2];

%

计算x值处的概率密度(cdf是累计分布)

y = pdf

(pd,x)

y

=

0.0540 0.2420 0.3989 0.2420 0.0540

同样，另一种表达

y = pdf(pd,x)

y =

0.0540 0.2420 0.3989 0.2420 0.0540

t分布概率密度函数tpdf(Student's

t probability density function)

语法

y =

tpdf(x,nu)

举例

tpdf(0,1:6)

ans =

0.3183 0.3536 0.3676 0.3750 0.3796 0.3827

相反，还可以通过p求t分布的t值

tinv (Student's t inverse cumulative distribution

function)

语法

x = tinv(p,nu)

举例

% the 99th percentile of the Student's t distribution for one to

six degrees of freedom

percentile = tinv(0.99,1:6)

percentile =

31.8205 6.9646 4.5407 3.7469 3.3649 3.1427

有一个问题，Matlab有一个inv矩阵求逆函数，不知与tinv什么关系，莫非tinv是在t分布下调用了inv计算程序？但p并不等是t的逆矩阵啊(即t*p

= E)啊？求解答

inv是矩阵求逆的意思。具体用法A=inv(B)，其中B是输入的可逆矩阵，输出A就是B的逆矩阵，逆矩阵满足性质 AB=BA=E

(E是单位阵)。如果输入的是不可逆矩阵会弹出警告，并返回inf。

调用举例：

>> inv([1 0;0 0])

警告: 矩阵为奇异工作精度。

ans =

Inf Inf

Inf Inf

>> inv(rand(2))

ans =

-13.0929 5.2640

12.0501 -3.3159

另附官方英文解释(输入doc inv也可以自己查看)：

Y = inv(X) returns theinverse of the square matrix X. A warning

messageis printed if X is badly scaled or nearly singular.

In practice, it is seldom necessary to form the explicit

inverseof a matrix. A frequent misuse of inv arises whensolving the

system of linear equations Ax = b.One way to solve this is with x =

inv(A)*b.A better way, from both an execution time and numerical

accuracy standpoint,is to use the matrix division operator x =

Ab.This produces the solution using Gaussian elimination, without

formingthe inverse. See mldivide ()for further information.