Well, I don't feel like working out the details of the algorithm now, but it seems like what you are looking for is the adverb /\. .
For a function f, we have f/ y_0 , y_1 , y_2 , y_3 --> y_0 f (y_1 f (y_2 f y_3)) . This can be used in many ways as a for loop: if you have a fixed list to loop over, and a verb update which makes updates your current value y based on the left argument x, then (update/ (list of values),init) loops over x and applies each to the current value using update. /\. is the same except that it stores the partial results to form a result of the same length as the input, for example, +/\. 4 3 2 1 10 6 3 1 So in general for this type of problem you would make such a verb update that depends on x_i and updates (m_i,s_i). It looks like you would also want to reverse the input before and after using &.|. . There is one problem in using this approach in that each x_i is a scalar and (m_i,s_i) is two values and must be stored in a list. That would be fixed if there were a dyadic form of / so that x_0 x_1 x_2 f y --> x_0 f x_1 f x_2 f y , and I have been missing such a verb many times. There are two equally ugly solutions for this: 1. Boxing: > f&.>/\. (<"_1 y) , <(m_init,s_init) 2. Making your own version: reduce = 1 : 'for_z. |.x do. y=.z u y end. y' NB. untested f reduce&(m_init,s_init) \. y The second case here doesn't allow for the optimized performance of /\. , and actually computes f reduce&(m_init,s_init) for each prefix of y, giving it time O(n^2). Therefore the first is the only viable option unless you want to write an adverb reduce_prefix which replicates /\. with a specified initial value. Marshall On Tue, Sep 13, 2011 at 5:47 PM, Johann Hibschman <[email protected]>wrote: > That's a different algorithm, one with poor numeric properties. For > example, > > %:((*<:)@#%~(#*+/@:*:)-*:@:(+/)) 1e6+?1e6#0 > 0.413838 > %:((*<:)@#%~(#*+/@:*:)-*:@:(+/)) 1e6+?1e6#0 > 0.145163 > > compared with the stddev from stats, which does a direct sum of squared > deviance: > > stddev 1e6+?1e6#0 > 0.288728 > stddev 1e6+?1e6#0 > 0.28884 > > Welford's method should have better properties than the sum of squared > deviance, but it has a strong sequential-update flavor that I don't know > how to do in J. > > P.S. As should be clear, I dropped a square root from my stddev > calculation below. > > "R.E. Boss" <[email protected]> writes: > > > ((*<:)@#%~(#*+/@:*:)-*:@:(+/)) > > > > > > R.E. Boss > > > > > >> -----Oorspronkelijk bericht----- > >> Van: [email protected] [mailto:programming- > >> [email protected]] Namens Johann Hibschman > >> Verzonden: dinsdag 13 september 2011 21:23 > >> Aan: [email protected] > >> Onderwerp: [Jprogramming] Welford's method for standard deviations > >> > >> I just came across Welford's method for calculating standard deviations, > >> and I realized I didn't know how to implement it in J. > >> > >> See > >> > >> - http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance > >> - http://www.johndcook.com/blog/2008/09/26/comparing-three-methods- > >> of-computing-standard-deviation/ > >> - http://www.johndcook.com/standard_deviation.html > >> > >> for the various articles that got me wondering about this. > >> > >> Here's the algorithm, as stated in the above links. > >> > >> 1. Initialize m_1 = x_1 and s_1 = 0. > >> 2. For subsequent xs, > >> m_k = m_{k-1} + (x_k - m_{k-1})/k > >> s_k = s_{k-1} + (x_k - m_{k-1})*(x_k - m_k) > >> 3. Then the kth running standard deviation is: > >> stddev_k = s_k / (k - 1) > >> > >> I've not tried to switch it from 1-based indexing to 0-based indexing, > >> or to use J-like phrasing in the problem, for fear I'd make a stupid > >> mistake. > >> > >> How would I do this in J? > >> > >> Cheers, > >> Johann > >> > >> ---------------------------------------------------------------------- > >> For information about J forums see http://www.jsoftware.com/forums.htm > > > > ---------------------------------------------------------------------- > > For information about J forums see http://www.jsoftware.com/forums.htm > > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
