require 'plot numeric png' NB. Gamma function and Beta function gamma=:!&<:! beta=: ] %@* [ !&<: + NB. J Vocabraly NB. Hypergeometrics NB. F is Fooks F=: 3 : 0 NB. Convert F(a;b;c;z) into monadic H. call 'a b c'=. 3{.y z=. > 3}.y m=. a,b n=. c m H. n z ) NB. Usage: F 1r2 ;1;3r2;z^2 NB. 1.01366 1.05912 1.15525 1.37327 NB. m=a+b : n=c NB. m H. n z=: }. 5%~ i.5 NB. 0.2 0.4 0.6 0.8 NB. Nishikawa Japla 1016/09 NB. Hypergeometric Function ================================== NB. Gauss Hypergeometric Function NB. test natural logarithm calc. NB. (1, 1, 2) Fn _0.5 for F(1, 1, 2;_0.5) NB. 0.810869 NB. NB. (^. (1 + 0.5)) % 0.5 for ln(1 + 0.5) % (0.5) NB. 0.81093 NB. (1, 1, 2) Fn 0.5 for F(1, 1, 2;0.5) NB. 1.38613 NB. (^. (1 - 0.5)) % (_0.5) for ln(1 - 0.5) % (_0.5) NB. 1.38629 Fn =: 3 : 0 : n =. 10 'a b c' =. x x0 =. y i =. 1 A =. a B =. b C =. c s =. 1 while. i < n do. NB. calc. coeficient ============================== NB. wr (2": i), ':', (3":A), '<', (":*/A), '> *', (3":B), '<', (":*/B),'> /',(3":C),'<', (":*/C),'>' CC =. ((,*/A) * (,*/B)) % (,*/C) NB. wr 'CC=', (": CC), ' : ', (": x: CC) NB. wr 's=', ":s A =. A, (a+i) B =. B, (b+i) C =. C, (c+i) NB. calc. x-value ================================= Xi =. (% ! i) * (x0^i) sx =. CC * Xi s =. s + sx i =. i + 1 end. s ) NB. Kummer Cofluent Hypergeometric Function ============ NB. (1, 1.5) Mn 0.5 NB. 1.41069 NB. (1, 1.5) Mn 1 NB. 2.03008 Mn =: 3 : 0 : n =. 10 'a c' =. x x0 =. y i =. 1 A =. a B =. 1 C =. c s =. 1 while. i < n do. NB. calc. coeficient NB. wr (2": i), ':', (3":A), '<', (":*/A), '> *', (3":B), '<', (":*/B),'> /',(3":C),'<', (":*/C),'>' CC =. ((,*/A) * (,*/B)) % (,*/C) NB. wr 'CC=', (": CC), ' : ', (": x: CC) NB. wr 's=', ":s A =. A, (a+i) B =. B, 1 NB. for Kummer Function NB. B =. B, (b+i) C =. C, (c+i) NB. calc. x-value Xi =. (% ! i) * (x0^i) sx =. CC * Xi s =. s + sx i =. i + 1 end. s ) NB. error function ============================ adj_erf =: 2p_0.5&* % ^ @: *: NB. J original =================== f1 =: (1 H. 1.5)@*: ferf =: f1 * adj_erf NB. ferf 0.5 1 NB. 0.5205 0.842701 NB. Nishikawa version 0 ============ m0 =: 3 : '(1, 1.5) Mn (*: y)' merf =: m0 * adj_erf NB. , merf"(0) 0.5 1 NB. 0.5205 0.842701 NB. Nishikawa version 1 ============ Merf =: 3 : 0 M1 =. (1 1.5) Mn *: y M2 =. ((2 % %: 1p1) * y) % (^ *: y) M =. M1 * M2 ) NB. Cumulative Normal Function 0-1 ============= CumNormf =: 3 : 0 0.5 * 1 + Merf (%: 0.5) * y ) NB. M.shimura NB. D =: N((^.S%X)+T*r(+,-)-:*:v)%v*%:T NB. M =: S,X*^-r*T BS =: monad define 'S X T r v' =. y -/(S,X*^-r*T)*N((^.S%X)+T*r(+,-)-:*:v)%v*%:T ) erf=:(*&(%:4p_1) % ^@:*:) * [: 1 H. 1.5 *: NB. Ewart Shaw Vector 18.4 cnd =:N=: [: -: 1: + [: erf %&(%:2) NB. CDF of N(0,1) a =: 60 65 0.25 0.08 NB. Probit NB. original is erf_rrobit.ijs NB. Japla 2010 require 'plot numeric trig files csv' require jpath '~addons/stats/distribs/normal.ijs' NB. ----------DATA ----------------- NB. Gujarati Basic Econometrics 4th ed. NB. (table 15.5) income(Xi, $1000) Ni housholders(ni) D0=: 6 40 8,8 50 12,10 60 18,:13 80 28 D0=: D0,15 100 45,20 70 36,25 65 39,:30 50 33 D0=: D0,35 40 30,:40 25 20 find_rate=: (%/)@(2 1 &{) "1 probit=:(%:@2:) * erfinv_pdistribs_@<:@+: NB. probit@find_rate D0 reg_probit=: probit@find_rate %. 1&,.@{."1 estim_probit=: 3 : '(reg_probit y)&p. {."1 y' probit_main=: pnormh_pdistribs_@estim_probit