---Devon McCormick wrote:
> Ambrus's solution is probably better than mine, but I did
> this all in J this
> way:
>
> NB.* inverseCND: inverse of cumulative norm. dist.; OK between 1e_10 &
> 1-1e_10.
> inverseCND=: 3 : 0"0
> if. (y>0.4999999998)*.y<0.500000002 do. 0 NB. Unstable
> near inflection
> else.
> f=: 13 : ((19j17":y),'-cnd y') NB. Function to solve
> f newton^:18]estInitInvCND y
> end.
> NB.EG _2.32634787404084-:inverseCND 0.01 NB. 1% of distribution to the
> left
> NB.EG 1.64485362695147-:inverseCND 0.95 NB. 5% of distribution to the right
> )
>
> ESTINVCND=: (i:6),:~1-cnd i:6 NB. Use to estimate inverse of cnd
>
> estInitInvCND=: 3 : 0"0
> whst=. 0 i.~/:y,0{ESTINVCND
> ix=. (<:1{$ESTINVCND)<.0>.whst-1 0
> -:+/ix{1{ESTINVCND
> )
>
> NB.* newton: Newton's iterative method to find zero of fnc
> given estimate.
> newton=: 1 : '- u % u D. 1'
>
Thanks Devon, I'll study this more closely. I also came up with a (simple)
iterative solution though not as succinct as yours. It also handles a mix of
sub-populations with different means & sds.
require 'stats/base/distribution'
NB.*pnorm01 v Probability of std normal above trunc point y
pnorm01=: [: -. normalcdf
NB.*zed01 v Height of std normal curve at truncation point y
zed01=: (%%:2p1) * ^@:(_0.5 * *:)
NB.*getTrunc v Trunc pt for top x members of sub-popln(s) y
NB. eg: _2.326347874-: 0.99 getTrunc 1 0 1 NB. 1% of distribution to the left
NB. eg: 1.644853627-: 0.05 getTrunc 1 0 1 NB. 5% of distribution to the right
NB. x is: number to select
NB. y is: 3-item list or 3-column array of sub-poplns to select from
NB. 0{"1 number in sub-poplns
NB. 1{"1 mean of sub-poplns
NB. 2{"1 stdev of sub-poplns
getTrunc=: 4 : 0
nsel=.x
'n mu sd'=. |:y
t=. {.mu
if. nsel>: +/n do. NB. select more than available
t-10*{.sd return.
end.
alpha=. 1
oldt=. t
oldnum=. nsel-~ +/ n * pnorm01 sd %~ t-mu NB. how many more needed
oldden=. - +/ n * sd %~ zed01 sd %~ t-mu
whilst. 1e_5 < | oldt-t do.
num=. nsel-~ +/ n * pnorm01 sd %~ t-mu
den=. - +/ n * sd %~ zed01 sd %~ t-mu
if. 0 > num * oldnum do. NB. too many selected
alpha=. -:alpha
t=. oldt - alpha*oldnum%oldden
else. NB. need to select more
'oldt oldnum oldden'=. t;num;den
alpha=. 1
t=. oldt - alpha*num%den
end.
end.
t
)
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
