Computing confidence intervals

Hi all.

I have this application where I have several minutes of
acoustic recordings at a few tens of kHz sampling frequency.
I need to compute power spectra from these data with
confidence intervals.

I use a one-pass variance estimator (pseudo code):

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X1=3D 0;             % Nfft x 1 data vectors
X2=3D 0;
for k=3D1:K
  Xk=3Dabs(fft(xk)); % Magnitude spectrum of k'th data frame
  X1=3DX1+xk;        % Accumulate magnitude spectra
  X2=3DX2+xk.^2;     % Accumulate power spectra
end

Xmean =3D X1/K;      % Compute the mean magnitude spectrum
Xvar  =3D X2/(K-1) - Xmean; % Compute variance of magnitude
                   % spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Now, I would like to compute various confidence intervals
based on the variance of the estimates. In order to do that,
I need the ERF function, defined implicitly through an
integral as something like

 x =3D ERF(p) =3D {x: 1/sqrt(2pi) \=ECnt_0^p exp(x^2)dx =3D p}

(I hope I got that one right...) and also the Q_nu(x) function
associated with the Chi-square distribution and Student's
t distribution.

The Q function is needed when one works with estimated
variances, as opposed to known variances.

I can hack the ERF function from available matlab commands,
but I can't see that basic matlab (no toolboxes) provides
a function that can be used to compute the Q function.

Does anyone know how to compute the Q function? I am tempted
to have a go at implementing it from scratch...

As I will need to implement the same application with C++
in a not too distant future, does anyone know the recipe for
the ERF function?

Any hints and suggestions are highly appreciated.

Rune

Rune wrote:

> (I hope I got that one right...) and also the Q_nu(x) function
> associated with the Chi-square distribution and Student's
> t distribution.
>
> The Q function is needed when one works with estimated
> variances, as opposed to known variances.
>
> I can hack the ERF function from available matlab commands,
> but I can't see that basic matlab (no toolboxes) provides
> a function that can be used to compute the Q function.
>
> Does anyone know how to compute the Q function? I am tempted
> to have a go at implementing it from scratch...
>
> As I will need to implement the same application with C++
> in a not too distant future, does anyone know the recipe for
> the ERF function?

You should be able to find all of that in numerical recipes, chapter 6:

http://www.library.cornell.edu/nr/cbookcpdf.html

There is a C++ version available, but not for download on the web. You
should be able to get it from a library, it's probably worth the
hassle, seeing that the C-version is quite old already.

Regards,
Andor

In article <1121343740.869828.104290@z14g2000cwz.googlegroups.com>,
 "Rune Allnor" <allnor@tele.ntnu.no> wrote:

> Hi all.
> 
> I have this application where I have several minutes of
> acoustic recordings at a few tens of kHz sampling frequency.
> I need to compute power spectra from these data with
> confidence intervals.
> 
> I use a one-pass variance estimator (pseudo code):
> 
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> X1= 0;             % Nfft x 1 data vectors
> X2= 0;
> for k=1:K
>   Xk=abs(fft(xk)); % Magnitude spectrum of k'th data frame
>   X1=X1+xk;        % Accumulate magnitude spectra
>   X2=X2+xk.^2;     % Accumulate power spectra
> end
> 
> Xmean = X1/K;      % Compute the mean magnitude spectrum
> Xvar  = X2/(K-1) - Xmean; % Compute variance of magnitude
>                    % spectrum
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> 
> Now, I would like to compute various confidence intervals
> based on the variance of the estimates. In order to do that,
> I need the ERF function, defined implicitly through an
> integral as something like
> 
>  x = ERF(p) = {x: 1/sqrt(2pi) \&#4294967295;nt 0^p exp(x^2)dx = p}
> 
> (I hope I got that one right...) and also the Q nu(x) function
> associated with the Chi-square distribution and Student's
> t distribution.
> 
> The Q function is needed when one works with estimated
> variances, as opposed to known variances.
> 
> I can hack the ERF function from available matlab commands,
> but I can't see that basic matlab (no toolboxes) provides
> a function that can be used to compute the Q function.


Sometimes you just need to reach up and pick the low
hanging fruit...

help erf

It would seem to take very little in the way of hacking.



> Does anyone know how to compute the Q function? I am tempted
> to have a go at implementing it from scratch...
> 
> As I will need to implement the same application with C++
> in a not too distant future, does anyone know the recipe for
> the ERF function?
> 
> Any hints and suggestions are highly appreciated.


If you really prefer to roll your own, then the first
source is always Abramowitz and Stegun. 

If this is not enough, then find a copy of Hart. It
appears to be out of print, but it is an absolute
treasure if you can find a copy. I was surprised that
Amazon did not find it. They did the last time I looked.

J.F. Hart, Computer Approximations (Wiley, New York, 1968)

HTH,
John D'Errico

John D'Errico wrote:
> In article <1121343740.869828.104290@z14g2000cwz.googlegroups.com>,
>  "Rune Allnor" <allnor@tele.ntnu.no> wrote:
>
> > Hi all.
> >
> > I have this application where I have several minutes of
> > acoustic recordings at a few tens of kHz sampling frequency.
> > I need to compute power spectra from these data with
> > confidence intervals.
> >
> > I use a one-pass variance estimator (pseudo code):
> >
> > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> > X1=3D 0;             % Nfft x 1 data vectors
> > X2=3D 0;
> > for k=3D1:K
> >   Xk=3Dabs(fft(xk)); % Magnitude spectrum of k'th data frame
> >   X1=3DX1+xk;        % Accumulate magnitude spectra
> >   X2=3DX2+xk.^2;     % Accumulate power spectra
> > end
> >
> > Xmean =3D X1/K;      % Compute the mean magnitude spectrum
> > Xvar  =3D X2/(K-1) - Xmean; % Compute variance of magnitude
> >                    % spectrum
> > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> >
> > Now, I would like to compute various confidence intervals
> > based on the variance of the estimates. In order to do that,
> > I need the ERF function, defined implicitly through an
> > integral as something like
> >
> >  x =3D ERF(p) =3D {x: 1/sqrt(2pi) \=ECnt 0^p exp(x^2)dx =3D p}
> >
> > (I hope I got that one right...) and also the Q nu(x) function
> > associated with the Chi-square distribution and Student's
> > t distribution.
> >
> > The Q function is needed when one works with estimated
> > variances, as opposed to known variances.
> >
> > I can hack the ERF function from available matlab commands,
> > but I can't see that basic matlab (no toolboxes) provides
> > a function that can be used to compute the Q function.
>
>
> Sometimes you just need to reach up and pick the low
> hanging fruit...
>
> help erf
>
> It would seem to take very little in the way of hacking.

Yes, I'm aware of it. Still, some hacking needs to be done,
ref the remarks on the erf documentation page. One needs the
NORMCDF form for computing confidence intervals.

> > Does anyone know how to compute the Q function? I am tempted
> > to have a go at implementing it from scratch...
> >
> > As I will need to implement the same application with C++
> > in a not too distant future, does anyone know the recipe for
> > the ERF function?
> >
> > Any hints and suggestions are highly appreciated.
>
>
> If you really prefer to roll your own, then the first
> source is always Abramowitz and Stegun.

Got that one...

> If this is not enough, then find a copy of Hart. It
> appears to be out of print, but it is an absolute
> treasure if you can find a copy. I was surprised that
> Amazon did not find it. They did the last time I looked.

I've seen it a few years ago. Seemed to be very useful.

> J.F. Hart, Computer Approximations (Wiley, New York, 1968)
>=20
> HTH,
> John D'Errico

Rune

"Rune Allnor" <allnor@tele.ntnu.no> wrote in message 
news:1121343740.869828.104290@z14g2000cwz.googlegroups.com...
>I need the ...... the Q_nu(x) function
>associated with the Chi-square distribution and Student's
>t distribution.

>The Q function is needed when one works with estimated
>variances, as opposed to known variances.

>Does anyone know how to compute the Q function? I am tempted
>to have a go at implementing it from scratch...



If by the Q function you mean the entity that is called Marcum's
Q function by radar folks, then there has been a lot of recent work
by M. K. Simon and M. S. Alouini on new representations and
calculation techniques for this function.  Much of this can be found
on IEEExplore if you have access to it.  Also, I believe they even have a
book on this subject.

Hope this helps.

Dilip V. Sarwate wrote:
> "Rune Allnor" <allnor@tele.ntnu.no> wrote in message
> news:1121343740.869828.104290@z14g2000cwz.googlegroups.com...
> >I need the ...... the Q_nu(x) function
> >associated with the Chi-square distribution and Student's
> >t distribution.
>
> >The Q function is needed when one works with estimated
> >variances, as opposed to known variances.
>
> >Does anyone know how to compute the Q function? I am tempted
> >to have a go at implementing it from scratch...
>
>
>
> If by the Q function you mean the entity that is called Marcum's
> Q function by radar folks,

I don't know. I am talking about the Q function that is tabulated in
most books on statistics to be used when computing confidence
intervals with estimated variance. I think the term "the chi-squared
distribution with nu degrees of freedom" is an often used term,
although I don't have a statistics textbook available to check the
details.

> then there has been a lot of recent work
> by M. K. Simon and M. S. Alouini on new representations and
> calculation techniques for this function.  Much of this can be found
> on IEEExplore if you have access to it.  Also, I believe they even have a
> book on this subject.
>
> Hope this helps.

It did. Lots of hits on IEEExplore. Will be useful, if we are talking
about the same thing.

Rune

"Rune Allnor" <allnor@tele.ntnu.no> wrote in message 
news:1121343740.869828.104290@z14g2000cwz.googlegroups.com...

Does anyone know how to compute the Q function? I am tempted
to have a go at implementing it from scratch...

As I will need to implement the same application with C++
in a not too distant future, does anyone know the recipe for
the ERF function?


Hello Rune,

I wrote some C code back in 1990 that calculates confidence intervals for 
Gaussian, Chi Square, Student's T and F Distributions. Would this be helpful 
to you? I took me a little while to find it, but I did and I've compiled it 
and it runs.

Clay

Clay S. Turner wrote:
> "Rune Allnor" <allnor@tele.ntnu.no> wrote in message
> news:1121343740.869828.104290@z14g2000cwz.googlegroups.com...
>
> Does anyone know how to compute the Q function? I am tempted
> to have a go at implementing it from scratch...
>
> As I will need to implement the same application with C++
> in a not too distant future, does anyone know the recipe for
> the ERF function?
>
>
> Hello Rune,
>
> I wrote some C code back in 1990 that calculates confidence intervals for
> Gaussian, Chi Square, Student's T and F Distributions. Would this be helpful
> to you? I took me a little while to find it, but I did and I've compiled it
> and it runs.
>
> Clay

That's exacly what I want to implement, but I would prefer to code it
myself for various legal reasons, as well as for the learning
experience.
I would be most interested in any references to the algorithm, though.

BTW, Clay, I get the impression you have coded a lot of neat stuff
at one time or another. Both the knowledge of how to do these things,
as well as the routines themselves, are the kind of "intellectual
capital"
that becomes more and more valuable as time goes by, since no one
does that kind of stuff anymore.

Perhaps something to keep in mind?

Rune

"Rune Allnor" <allnor@tele.ntnu.no> wrote in message 
news:1121515426.925225.45080@g49g2000cwa.googlegroups.com...
>
>
> Clay S. Turner wrote:
>> "Rune Allnor" <allnor@tele.ntnu.no> wrote in message
>> news:1121343740.869828.104290@z14g2000cwz.googlegroups.com...
>>
>> Does anyone know how to compute the Q function? I am tempted
>> to have a go at implementing it from scratch...
>>
>> As I will need to implement the same application with C++
>> in a not too distant future, does anyone know the recipe for
>> the ERF function?
>>
>>
>> Hello Rune,
>>
>> I wrote some C code back in 1990 that calculates confidence intervals for
>> Gaussian, Chi Square, Student's T and F Distributions. Would this be 
>> helpful
>> to you? I took me a little while to find it, but I did and I've compiled 
>> it
>> and it runs.
>>
>> Clay
>
> That's exacly what I want to implement, but I would prefer to code it
> myself for various legal reasons, as well as for the learning
> experience.
> I would be most interested in any references to the algorithm, though.
>
> BTW, Clay, I get the impression you have coded a lot of neat stuff
> at one time or another. Both the knowledge of how to do these things,
> as well as the routines themselves, are the kind of "intellectual
> capital"
> that becomes more and more valuable as time goes by, since no one
> does that kind of stuff anymore.
>
> Perhaps something to keep in mind?
>
> Rune
>

Hello Rune,

Since I wrote the code for myself - there are no legal issues on this end. I 
recall I had finished taking my stats courses and I had a project where I 
needed the intervals. Certainly looking in a table in the back of a stats 
book would have sufficed for the project, but I thought it would be neat to 
figure out how to do this. Basically I used Gauss-LeGendre quadrature 
integration coupled with Newton's iteration to find the intervals.

So imagine you have a probability density, f(x),  and you wish to find which 
"b" makes the integral from a to b of f(x)dx = p. "a" is usually equal to 
zero.  Then you just guess a value for b and use Newton's iteration as 
follows:

b_hat = b + (p - integral from a to b of f(x) dx ) /  f(b)

then use b_hat as b for the next iteration. This will only require a few 
steps to converge.

I used a 3pt Gauss-LeGendre quadrature integration formula to evaluate the 
above integral. This has the advantage of using only interior points and is 
a lot more accurate than Simpson's rule for basically the same effort.

I find that I often get called to solve problems that have arcane answers 
and having this knowledge is starting to land me work. I have written lots 
of code snippets over the years to figure out various problems. Often when 
I'm presented with a problem - I ask myself how would I program a computer 
to solve it and then I write it. I guess this is how I think about things.

Currently I'm working on my dissertation and most of my algorithmic effort 
is spent there these days. I'm developing a new technique for solving 
certain problems in atomic physics. The resulting code so far is a couple of 
orders of magnitude faster and much more accurate than other versions. This 
general algo will be applicable to molecular physics as well. For example 
when solving the atomic Hydrogen eigenvalue problem (of which the 
theoretical answers are known and may be used as a reference) my approach 
finds the energies correct to 14 decimal places using math that supports 
15.5 decimal places. I like working with numerical algorithms because so 
much is still to be discovered. Floating point noise presents interesting 
challenges.

But one day I may put together a lot of my code snippets into a tome. It 
will be eclectic though.

Clay

 p.s. Let me know if you have any more questions about how to find 
confidence intervals. If you recall most statisticians work with right hand 
tail probabilities. So you may want to complement the probability because 
you don't want to integrate out to infinity.

Rune wrote:

> Perhaps something to keep in mind?

Are you hiring?

:-)