[R] function of probability for normal distribution

[Greg Snow]

It looks like you are reinventing wheels here (not necessarily a bad thing).

If you want the probabilities associated with the normal distribution, then look at the help for dnorm and pnorm functions.

If you are recreating this as a learning experience/homework, then you should be up front about that, what is the assignment? What level of help are you allowed? Etc.  Then be more explicit in what parts you understand and what you need help with.  This list is generally not for homework help (definitely not for doing it for you), but there are cases where the question and limitations are clear that we have given appropriate hints.

-- 
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.