;+
; NAME:
;    multigauss
;
; PURPOSE:
;    This function estimates the mean and standard deviation of an ensemble of 
;    Gaussian probability density functions.
;
; CATEGORY:
;    Statistics
;
; CALLING SEQUENCE:
;    result = multigauss( in_mean, in_std )
;
; INPUTS:
;    IN_MEAN:  A required float vector of length N_DATA containing the means 
;        of the N_DATA Gaussian probability density functions with the standard 
;        deviations specified in IN_STD.
;    IN_STD:  A required float vector of length N_DATA containing the standard 
;        deviations of the N_DATA Gaussian probability density functions with 
;        the means specified in IN_MEAN.
;    N_LOC, LOC_RANGE
;
; KEYWORD PARAMETERS:
;    LOC_RANGE:  An optional 2-element float vector specifying the range of the 
;        N_LOC points in OUT_LOC at which to sample the ensemble probability 
;        density function.  The default is the the minimum/maximum of the mean 
;        minus/plus three times the standard deviation of the input density 
;        functions.
;    N_LOC:  An optional integer specifying the number of points at which to 
;        to sample the ensemble probability density function.  The default is 
;        1000.
;    OUT_CDF:  Returns a N_LOC float vector of the cumulative distribution 
;        function of the ensemble of Gaussian density functions estimated at 
;        the locations in OUT_LOC.
;    OUT_LOC:  Returns a N_LOC float vector specifying the locations at which 
;        to estimate the probability density function of the ensemble of 
;        Gaussian density functions, output as OUT_PDF, as well as the 
;        cumulative distribution function of the ensemble, output as OUT_CDF.
;    OUT_PDF:  Returns a N_LOC float vector of the probability density function 
;        of the ensemble of Gaussian density functions estimated at the 
;        locations in OUT_LOC.
;
; OUTPUTS:
;    RESULT:  A float vector of length 2 returning the mean and standard 
;        deviation of the ensemble of density functions.
;    OUT_CDF, OUT_LOC, OUT_PDF
;
; USES:
;    -
;
; PROCEDURE:
;    This function adds up the input probability density functions and then 
;    calculates the percentiles of the result at the Gaussian mean and 
;    plus/minus one standard deviation.
;
; EXAMPLE:
;    ; Estimate mean and standard deviation of the addition of two units 
;    ; Gaussian distributions, one centred at zero and the other at one.
;    print, multigauss( [0.,1.], [1.,1.] )
;    ; The result should be about [ 0.5, 1.125 ].
;
; MODIFICATION HISTORY:
;    Written by:  Daithi A. Stone (dastone@runbox.com), 2026-06-24
;-

;***********************************************************************

FUNCTION MULTIGAUSS, $
    IN_MEAN, IN_STD, $
    N_LOC=n_loc, $
    OUT_CDF=out_cdf, OUT_LOC=out_loc, OUT_PDF=out_pdf

;***********************************************************************
; Constants and options

; The number of input data points
n_in = n_elements( in_mean )
if n_in eq 0 then stop
if n_elements( in_std ) ne n_in then stop

; The default range of values to sample
if not( keyword_set( loc_range ) ) then begin
  loc_range = [ min( in_mean - 3 * in_std ), max( in_mean + 3 * in_std ) ]
endif
; The default number of sampling points
if not( keyword_set( n_loc ) ) then n_loc = 1000

;***********************************************************************
; Estimate parameters

; Generate vector of points to sample
out_loc = loc_range[0] $
    + findgen( n_loc ) / ( n_loc - 1. ) * ( loc_range[1] - loc_range[0] )
; Initialise ensemble density function vector
out_pdf = fltarr( n_loc )

; Iterate through Gaussian functions
for i_in = 0, n_in - 1 do begin
  ; Generate this probability density function
  temp_pdf = exp( -0.5 * ( ( out_loc - in_mean[i_in] ) / in_std[i_in] ) ^ 2. )
  temp_pdf = temp_pdf / total( temp_pdf )
  ; Add this function to the ensemble
  out_pdf = out_pdf + temp_pdf
endfor
; Normalise the ensemble function
out_pdf = out_pdf / total( out_pdf )
; Convert to a cumulative distribution function
; (Subtracting the first value avoids the low-bias issue when normalising by 
; the maximum value.)
out_cdf = total( out_pdf, cumulative=1 ) - out_pdf[0]
out_cdf = out_cdf / max( out_cdf )

; Calculate the mean and standard deviation of this density function
temp = interpol( out_loc, out_cdf, gauss_pdf([-1.,0.,1.]) )
params = [ temp[1], ( temp[2] - temp[0] ) / 2. ]

;***********************************************************************
; The end

return, params
END