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