[R] function of probability for normal distribution

[Rene]

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.