from org.virbo.dataset.DataSetOps import moment
from math import sqrt

len= 100000

tmu= 2.24
tsd= 0.80

dist= randn( len ) * tsd + tmu
r= where( dist < -1 )
dm2= dist[r]
r= where( ( dist >= -1 ) & ( dist < 0 ) )
dm1= dist[r]
r= where( ( dist >=  0 ) & ( dist < 1 ) )
dp1= dist[r]
r= where( ( dist >=  1 ) )
dp2= dist[r]


mu= zeros(4)
sd= zeros(4)
vv= zeros(4)
nn= zeros(4)

nn[0]= dm2.length()
nn[1]= dm1.length()
nn[2]= dp1.length()
nn[3]= dp2.length()

mu[0]= moment(dm2).value()
mu[1]= moment(dm1).value()
mu[2]= moment(dp1).value()
mu[3]= moment(dp2).value()

sd[0]= moment(dm2).property( 'stddev' )
sd[1]= moment(dm1).property( 'stddev' )
sd[2]= moment(dp1).property( 'stddev' )
sd[3]= moment(dp2).property( 'stddev' )

vv= sd**2

m= moment(dist)
mean= m.value()

# correct each of the four populations to the actual mean
vv= ( nn*vv + nn * ( mu - mean )**2 ) / nn

onepassvar= total( vv*nn ) / dist.length()
onepassmean= total( mu*nn ) / dist.length()

print 'true', tmu, tsd**2, tsd
print 'twopass', mean, m.property( 'stddev' )**2, m.property( 'stddev' )
print 'onepass', onepassmean, onepassvar, sqrt(onepassvar)