-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathqnorm.c
67 lines (57 loc) · 1.3 KB
/
qnorm.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
/* Ratfor code of qnorms, from S:
real function qnorms(pr)
real pr
real p,eta,term,f1,f2,f3,f4,f5,f6
data f1/.010328/
data f2/.802853/
data f3/2.515517/
data f4/.001308/
data f5/.189269/
data f6/1.432788/
if(pr<=0.|pr>=1.)ERROR(`Probability out of range (0,1)')
p = pr
if(p > 0.5) p = 1.0 - pr
# depending on the size of pr this may error in alog or sqrt
eta=sqrt(-2.*alog(p))
term=((f1*eta+f2)*eta+f3)/(((f4*eta+f5)*eta+f6)*eta+1.0)
if(pr<=.5)return( term - eta )
else return ( eta - term )
end
*/
/*
cc -o qnorm qnorm.c
*/
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
/* adapted from Ratfor code used in S */
double
qnorm(double pr)
{
double p, eta, term,
f1 = .010328,
f2 = .802853,
f3 = 2.515517,
f4 = .001308,
f5 = .189269,
f6 = 1.432788;
if(pr <= 0. || pr >= 1.) printf("Probability out of range (0,1): %f", pr);
p = pr;
if(p > 0.5) p = 1.0 - pr;
/* depending on the size of pr this may error in
log or sqrt */
eta = sqrt(-2.0*log(p));
term =((f1*eta+f2)*eta+f3)/(((f4*eta+f5)*eta+f6)*eta+1.0);
if(pr <= .5) return( term - eta );
else return ( eta - term );
}
main()
{
int i;
double p;
for(i=0; i<100; i++) {
printf("Enter a probability: ");
scanf("%lf", &p);
printf(" normal quantile = %f \n", qnorm(p));
}
}