Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.snow at imail.org
801.408.8111
> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-
> project.org] On Behalf Of Rene
> Sent: Thursday, August 20, 2009 5:19 AM
> To: r-help at r-project.org
> Subject: Re: [R] function of probability for normal distribution
>
> Dear all,
>
> Because my last email was in html format, so it was a disaster to read.
> I
> have second thought of my question asked in my last email, and came up
> some
> solution to myself, but I found the result was a bit weird, can someone
> please help look at my coding and advise where I have done wrong?
>
> I need to write a function which compute the probability for standard
> normal
> distribution. The formula is
>
> P(z)= 1/2 + 1/(sqrt(2*pi)) * sum ((-1)^k*z^(2*k+1) /(2^k*k!*(2*k+1))
>
>
> >From series expansion for the exponential function, we know e^x= sum
> (x^n/n!), we can substitute x = -t^2/2, so e^(-t^2/2) = sum
> ((-1)^k*t^(2*k)/(2^k*k!)). Comparing with the formula above, the sum
> ((-1)^k*z^(2*k+1) /(2^k*k!*(2*k+1)) is very similar to e^(-t^2/2),
> except
> there is z/(2*k+1) in the formula.
>
> So I create a function for the sum ((-1)^k*z^(2*k+1) /(2^k*k!*(2*k+1)),
> which is
>
> expf = function (t)
> {
> neg=(t<0)
> t = ifelse(neg,-t,t)
> k=0;term=1;sum=1
> oldsum=0
> while(any(sum!=oldsum)){
> oldsum = sum
> k = k+1
> term = term*(((-1)^1*t^2) / (2*k))
> sum = sum + (t/(2*k+1))*term
> }
> ifelse(neg,1/sum,sum)
> }
>
> Now the sum part of the probability can be replaced by expf function,
> then I
> create the function for the probability p(z), which is
>
> prob = function(z)
> {
> 1/2+(1/sqrt(2*pi))*expf(z)
> }
>
> Am I doing the right thing here? However, when I test the probability
> function, e.g. prob(1:10), the result appear some negative value and
> some
> very large value which do not seem right to me as probability values.
> Can
> someone please guide me where I have done wrong here?
>
>
> Thanks a lot.
>
> Rene.
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-
> guide.html
> and provide commented, minimal, self-contained, reproducible code.