Rene
2009-Aug-19 11:04 UTC
[R] R function for Probabilities for the (standard) normal distribution
Dear All, I need to write an R function which computes values of Probabilities for the (standard) normal distribution, ¦µ(z) for z > 0 by summing this power series. (We should accumulate terms until the sum doesn't change). I also need to make the function work for vector of values z. The initial fomular is ¦µ(z) = ( 1/ sqrt(2¦Ð) )* ¡Ò e^(-t^2/2)dt (¡Ò is from -¡Þ, z) = 1/ 2 + ( 1/ sqrt(2¦Ð) )* ¡Ò e^(-t^2/2)dt (¡Ò is from 0 to z) I can substituting x = -t^2/2 into the series expansion for the exponential function e^x = ¡Æ x^n/n! (¡Æ is from n=0 to ¡Þ) I can obtain the series e^(-t^2/2) = ¡Æ (-1)^k*t^2k / 2^k*k! (¡Æ is from n=0 to ¡Þ) This series can be integrated term by term to obtain the formula ¦µ(z) = 1/ 2 + ( 1/ sqrt(2¦Ð) ) * ¡Æ (-1)^k*z^(2k+1) / (2^k*k! *(2k+1)) (¡Æ is from n=0 to ¡Þ) I know how to write the R function for exponential function e^x expf = function (x) { x=ifelse((neg=(x<0)),-x,x) n=0;term=1 oldsum=0; newsum=1 while(any(newsum != oldsum)) { oldsum=newsum n=n+1 term = term*x/n newsum = newsum+term} ifelse(neg, 1/newsum, newsum) } I know it will be similar to the above coding, but I don¡¯t know exactly how should we modify the above coding in order to get Probabilities for the (standard) normal distribution, ¦µ(z) for z > 0. Can anybody advise me on this?? Thanks a lot. Rene. [[alternative HTML version deleted]]
David Winsemius
2009-Aug-19 11:38 UTC
[R] R function for Probabilities for the (standard) normal distribution
On Aug 19, 2009, at 7:04 AM, Rene wrote:> Dear All, > > I need to write an R function which computes values of > Probabilities for > the (standard) normal distribution, ??(z) for z > 0 by summing this > power > series. (We should accumulate terms until the sum doesn't change). I > also > need to make the function work for vector of values z.Two problems: One, this appears to be a homework exercise.> > The initial fomular is > > ??(z) = ( 1/ sqrt(2??) )* ?? e^(-t^2/2)dt (?? is from -??, z) > > = 1/ 2 + ( 1/ sqrt(2??) )* ?? e^(-t^2/2)dt (?? is from 0 > to z)Two, you are posting in HTML mail and the characters that you may think we should be seeing reading as Greek letters including I presume capital sigma for summation are .... not displayed properly.> > > > I can substituting x = -t^2/2 into the series expansion for the > exponential > function > > e^x = ?? x^n/n! (?? is from n=0 to ??) > > I can obtain the series > > e^(-t^2/2) = ?? (-1)^k*t^2k / 2^k*k! (?? is from n=0 to ??) > > This series can be integrated term by term to obtain the formula > > ??(z) = 1/ 2 + ( 1/ sqrt(2??) ) * ?? (-1)^k*z^(2k+1) / (2^k*k! *(2k > +1)) > (?? is from n=0 to ??) > > > I know how to write the R function for exponential function e^x > > expf = function (x) > > { > > x=ifelse((neg=(x<0)),-x,x) > > n=0;term=1 > > oldsum=0; newsum=1 > > while(any(newsum != oldsum)) { > > oldsum=newsum > > n=n+1 > > term = term*x/n > > newsum = newsum+term} > > ifelse(neg, 1/newsum, newsum) > > } > > I know it will be similar to the above coding, but I don??t know > exactly how > should we modify the above coding in order to get Probabilities for > the > (standard) normal distribution, ??(z) for z > 0. > > Can anybody advise me on this?? > > Thanks a lot. > > Rene. > [[alternative HTML version deleted]] > > ______________________________________________David Winsemius, MD Heritage Laboratories West Hartford, CT
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.
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 > stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide R-project.org/posting- > guide.html > and provide commented, minimal, self-contained, reproducible code.