作者markoo (上行下笑阿..)
看板NCTU-STAT98G
標題[情報] 參考資料
時間Thu Aug 13 13:28:39 2009
#define PI 3.1415926
#define f(x) exp(-x*x/2)/sqrt(2*PI)
#define N 1000
double integral(double z, int size)
{
int i;
double up,low,avg,ans=0;
for (i=0;i<size;i++)
{
up=f((i*z*1.0)/size);
low=f(((i+1)*z*1.0)/size);
avg=(up+low)*(z/size)/2;
ans+=avg;
}
return ans;
}
double cdf(double z)//usage: P{Z<=z} = p, Z~N(0,1), input z.
{
double ans=0;
if (z<0) ans=0.5-integral(-z,N);
else ans=0.5+integral(z,N);
return ans;
}
double icdf(double prob,double precise)
//usage: P{Z<=z} = p, Z~N(0,1), input p and precise (精確度:建議0.001即可)
{
double ans=0;
if (prob>=0.5)
{
prob=prob-0.5;
do {
ans+=precise;
} while (integral(ans,N) <=prob);
}
else
{
prob = 0.5-prob;
do {
ans+=precise;
} while (integral(ans,N)<=prob);
ans = -ans;
}
return ans;
}
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 140.113.7.248
※ 編輯: markoo 來自: 140.113.7.248 (08/13 13:38)