R help - Jun 2006 - Standard Deviation Distribution

I'm having trouble with the standard deviation distribution
as shown on http://mathworld.wolfram.com/StandardDeviationDistribution.html .
(Eric Weisstein references Kenney and Keeping 1951, which I can't check.)

I believe the graphs they show, but when I code the function in R, according to
the listed formula,
I get very different graphs.

Would someone please point out my error or tell me where it's already
implemented in R?

Here is my version:
> sddist
function(s,n) {
sig2 <- n*s*s/(n-1)
2*(n/(2*sig2))^((n-1)/2) / gamma((n-1)/2) * exp(-n*s*s/(2*sig2)) * s^(n-2)
}

Version 2.3.1 (2006-06-01) on Windows XP SP2

Thanks for any help.

David L. Reiner
Rho Trading Securities, LLC
Chicago? IL? 60605
312-362-4963

davidr at rhotrading.com wrote:
> I'm having trouble with the standard deviation distribution
> as shown on http://mathworld.wolfram.com/StandardDeviationDistribution.html
.
> (Eric Weisstein references Kenney and Keeping 1951, which I can't
check.)
>
> I believe the graphs they show, but when I code the function in R,
according to the listed formula,
> I get very different graphs.
>
> Would someone please point out my error or tell me where it's already
implemented in R?
>
> Here is my version:
>   
>> sddist
>>     
> function(s,n) {
> sig2 <- n*s*s/(n-1)
> 2*(n/(2*sig2))^((n-1)/2) / gamma((n-1)/2) * exp(-n*s*s/(2*sig2)) * s^(n-2)
> }
>
>   
Your initial factor looks a lot different from theirs.  I think you 
believed your variable name (i.e. sig2 is sigma^2), but didn't
define it that way.

Duncan Murdoch
> Version 2.3.1 (2006-06-01) on Windows XP SP2
>
> Thanks for any help.
>
> David L. Reiner
> Rho Trading Securities, LLC
> Chicago  IL  60605
> 312-362-4963
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
>

davidr at rhotrading.com wrote:
> > sddist
>>> > >>
> > function(s,n) {
> > sig2 <- n*s*s/(n-1)
> > 2*(n/(2*sig2))^((n-1)/2) / gamma((n-1)/2) * exp(-n*s*s/(2*sig2)) *
s^(n-2)
> > }

There is another way using more statistical knowledge:
Building on functions already available in R
one may note that the easiest way of defining
sddist is

sddist<-function(s,n) dchisq(n*s^2,n-1)*2*n*s

Plotting this function for n in seq(2,12,2)
reproduces the graph in Mathworld.

The idea is the following:
given a random variable x with a chisquare distribution with
df n-1, s can be defined by s=sqrt(x/n) and therefore
x=n*s^2

If F and G are the cumulative distribution functions of x and s
and f and g are the density functions of x and s, then we have
G(s)=F(n*s^2) and therefore
g(s)=f(n*s^2)*2*n*s
-- 
Erich Neuwirth, University of Vienna
Faculty of Computer Science
Computer Supported Didactics Working Group
Visit our SunSITE at http://sunsite.univie.ac.at
Phone: +43-1-4277-39464 Fax: +43-1-4277-39459