;====================================================================
;; icFamStatCore
;; collect non-central moments of 1st..4th order
;====================================================================
procedure( icFamStatCore(wave work)
	let( (xVec yVec len y1 y2 y3 y4 N v p s)

		xVec = drGetWaveformXVec(wave)
		yVec = drGetWaveformYVec(wave)
		len  = drVectorLength(xVec)

		if( !work->xVec then
			;; initialize
						
			y1= drCreateVec('double len)
			y2= drCreateVec('double len)
			y3= drCreateVec('double len)
			y4= drCreateVec('double len)
						
			for(i 0 len-1
				v = drGetElem(yVec i)

				p = v
				drAddElem(y1 p)
				p = p*v
				drAddElem(y2 p)
				p = p*v
				drAddElem(y3 p)
				p = p*v
				drAddElem(y4 p)
			)
			y1->units  = yVec->units
			y1->name   = "mean"
			y2->units  = yVec->units
			y2->name   = "sqr"
			y3->units  = yVec->units
			y3->name   = "cub"
			y4->units  = yVec->units
			y4->name   = "quad"

			work->xVec = xVec
			work->y1   = y1
			work->y2   = y2
			work->y3   = y3
			work->y4   = y4
			work->N    = makeVector(len 1)
			work->len  = len
		else
			when( xVec!=work->xVec
				error("icFamStat: different x-axis in %L\n" wave)
			)

			for(i 0 len-1
				v = drGetElem(yVec i)
				when(v 
					N = work->N[i] + 1
					work->N[i] = N

					p = v
					;; recurrence formulae (mean of y)
					s = drGetElem(work->y1 i)
					s = s + (p - s) / float(N)
					drSetElem(work->y1 i s)

					p = p*v
					;; recurrence formulae (mean of y^2)
					s = drGetElem(work->y2 i)
					s = s + (p - s) / float(N)
					drSetElem(work->y2 i s)

					p = p*v
					;; recurrence formulae (mean of y^3)
					s = drGetElem(work->y3 i)
					s = s + (p - s) / float(N)
					drSetElem(work->y3 i s)

					p = p*v
					;; recurrence formulae (mean of y^4)
					s = drGetElem(work->y4 i)
					s = s + (p - s) / float(N)
					drSetElem(work->y4 i s)
				)
			)
		); if
	); let
)
;====================================================================

;====================================================================
;; icFamStatHop
;; iterate over families, call core with single waveforms
;; collect data in work
;====================================================================
procedure( icFamStatHop(wave work)
	let( (wf)
		cond(
			( drIsWaveform(wave)
				 icFamStatCore(wave work)
			)
			( famIsFamily(wave)
				;; call self with all family members
				foreach(sv famGetSweepValues(wave)
					wf = famValue(wave sv)
					icFamStatHop(wf work)
				)
			)
			( t
				error("icFamStat: cannot handle %L\n" wave)
			)
		); cond
	); let
)
;====================================================================

;====================================================================
;; icFamStat(type wave(s))
;; do statistic over a collection of waves/families
;====================================================================
procedure( icFamStat(type @rest waves)
	let( (work wf fm m v s k N scale)
		work = ncons(nil)
		
		;; run over all input waveforms
		foreach(wave waves
			icFamStatHop(wave work)
		)

		;; build result from collected data
		if( type=="non-central moments" then
			;; EX + non-central moments EX^2, EX^3, EX^4
			fm = famCreateFamily("statistics" 'string)
			wf = drCreateWaveform(work->xVec work->y1)
			famAddValue(fm "1st" wf)
			wf = drCreateWaveform(work->xVec work->y2)
			famAddValue(fm "2nd" wf)
			wf = drCreateWaveform(work->xVec work->y3)
			famAddValue(fm "3rd" wf)
			wf = drCreateWaveform(work->xVec work->y4)
			famAddValue(fm "4th" wf)

		else
			;; central/standardized moments
			;; sigma^2 = variance = E(X^2-u^2) = E(X^2) - E(X)^2
			;; skew     = E( ((X-u)/sigma)^3 ) = ( E(X^3) - 3*u*sigma^2 - u^3 ) / sigma^3
			;; kurtosis = E( ((X-u)/sigma)^4 ) = (E(X^4) - 4*u*sigma^3*skew - 6*u^2*sigma^2 - u^4) / sigma^4

			for(i 0 work->len-1

				; for unbiased variance: rescale from 1/N to 1/N-1 
				; better estimate for N>6, otherwise use biased variance
				N = work->N[i]
				scale = if(N<6 then 1.0 else N / float(N-1) )
				
				m = drGetElem(work->y1 i)
				v = drGetElem(work->y2 i)
				s = drGetElem(work->y3 i)
				k = drGetElem(work->y4 i)

				v = (v - m*m) * scale 
				if( v<=0 then
					v = 0
					s = 0
					k = 0
				else
					s = s - 3*m*v - m*m*m
					k = k - 4*m*s - 6*m*m*v - m*m*m*m
					
					unless( type=="central moments"
						v = sqrt(v)
						s = s / (v*v*v)
						k = k / (v*v*v*v)
					)
				)

				drSetElem(work->y2 i v)
				drSetElem(work->y3 i s)
				drSetElem(work->y4 i k)
			); for
			
			if( type=="central moments" then
				work->y2->name  = "variance"
				work->y3->name  = "u3"
				work->y4->name  = "u4"
			else
				;; standardized moments
				work->y2->name  = "stdev"
				work->y3->name  = "skew"
				work->y4->name  = "kurtosis"
			)
			
			fm = famCreateFamily("statistics" 'string)
			wf = drCreateWaveform(work->xVec work->y1)
			famAddValue(fm work->y1->name wf)
			wf = drCreateWaveform(work->xVec work->y2)
			famAddValue(fm work->y2->name wf)
			wf = drCreateWaveform(work->xVec work->y3)
			famAddValue(fm work->y3->name wf)
			wf = drCreateWaveform(work->xVec work->y4)
			famAddValue(fm work->y4->name wf)
		); if
		fm
	); let
)
;====================================================================


;====================================================================
;; GUI builder information
;====================================================================
ocnmRegGUIBuilder(
 '(nil
  function icFamStat
  name icFamStat
  description "Statistics over family graphs"
  category ("Custom Functions")
  analysis (nil
	general (nil
		args (type wave)
		signals (nil
			wave (nil
				prompt "Signal"
				tooltip "Signal"
			)
		)
		params(nil
			type (nil
				prompt "type"
				tooltip "statistics"
				guiRowHint 1
				default "standardized moments"
				type (
					"standardized moments"
					"central moments"
					"non-central moments"
				)
				required nil
			)
		)
		inputrange t
	)
  )
  outputs(result)
 )
)
;====================================================================