{"action":"created","issue":{"url":"https://api.github.com/repos/stan-dev/stan/issues/1157","labels_url":"https://api.github.com/repos/stan-dev/stan/issues/1157/labels{/name}","comments_url":"https://api.github.com/repos/stan-dev/stan/issues/1157/comments","events_url":"https://api.github.com/repos/stan-dev/stan/issues/1157/events","html_url":"https://github.com/stan-dev/stan/issues/1157","id":51197449,"number":1157,"title":"Inverse CDF for normal","user":{"login":"rtrangucci","id":3925639,"avatar_url":"https://avatars.githubusercontent.com/u/3925639?v=3","gravatar_id":"","url":"https://api.github.com/users/rtrangucci","html_url":"https://github.com/rtrangucci","followers_url":"https://api.github.com/users/rtrangucci/followers","following_url":"https://api.github.com/users/rtrangucci/following{/other_user}","gists_url":"https://api.github.com/users/rtrangucci/gists{/gist_id}","starred_url":"https://api.github.com/users/rtrangucci/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/rtrangucci/subscriptions","organizations_url":"https://api.github.com/users/rtrangucci/orgs","repos_url":"https://api.github.com/users/rtrangucci/repos","events_url":"https://api.github.com/users/rtrangucci/events{/privacy}","received_events_url":"https://api.github.com/users/rtrangucci/received_events","type":"User","site_admin":false},"labels":[{"url":"https://api.github.com/repos/stan-dev/stan/labels/Feature","name":"Feature","color":"84b6eb"}],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":11,"created_at":"2014-12-06T21:51:47Z","updated_at":"2015-03-08T13:40:06Z","closed_at":null,"body":"Written by @bgoodri: Stan has the CDF and inverse CDF of the logistic distribution, but lacks the inverse CDF of the standard normal distribution (and many other distributions) in C++. Remember that the derivative of the inverse CDF is the reciprocal of the density (including constants).\r\n\r\nThere is a rational approximation algorithm for the inverse of the standard normal CDF here:\r\nhttp://home.online.no/~pjacklam/notes/invnorm/\r\nwhich can be implemented (and tested) in Stan code with\r\n```\r\n/** Stan program to demonstrate useage of the standard normal inverse CDF.\r\n* Algorithm first derived in 2003 by Peter Jon Aklam at \r\n* http://home.online.no/~pjacklam/notes/invnorm/\r\n* @param p Number(s) on the (0,1) interval\r\n* @param a Vector of constants defined in the transformed data block\r\n* @param b Vector of constants defined in the transformed data block\r\n* @param c Vector of constants defined in the transformed data block\r\n* @param d Vector of constants defined in the transformed data block\r\n* @param N Integer defining length of vector or array in non-scalar case\r\n* @return Number(s) on the entire real line, almost equal to Phi^{-1}(p)\r\n*/\r\n\r\nfunctions {\r\n  real invnorm(real p, vector a, vector b, vector c, vector d) {\r\n    real x; real q; real r; real p_low;  real p_high; real e; real u;\r\n    if( p <= 0 )     return(negative_infinity());\r\n    else if (p >= 1) return(positive_infinity());\r\n    p_low  <- 0.02425;\r\n    p_high <- 0.97575; # 1 - p_low\r\n    if ( (p_low <= p) && (p <= p_high) ) { # central region\r\n      q <- p - 0.5;\r\n      r <- square(q);\r\n      x <- (((((a[1]*r+a[2])*r+a[3])*r+a[4])*r+a[5])*r+a[6])*q /\r\n           (((((b[1]*r+b[2])*r+b[3])*r+b[4])*r+b[5])*r+1.0);      \r\n    }\r\n    else if (p < p_low) { # lower region\r\n      q <- sqrt(-2.0*log(p));\r\n      x <- (((((c[1]*q+c[2])*q+c[3])*q+c[4])*q+c[5])*q+c[6]) /\r\n            ((((d[1]*q+d[2])*q+d[3])*q+d[4])*q+1.0);            \r\n    }\r\n    else { # upper region\r\n      q <- sqrt(-2.0*log1m(p));\r\n      x <- -(((((c[1]*q+c[2])*q+c[3])*q+c[4])*q+c[5])*q+c[6]) /\r\n             ((((d[1]*q+d[2])*q+d[3])*q+d[4])*q+1.0);\r\n    }\r\n\r\n    # refinement\r\n    e <- Phi(x) - p;\r\n    u <- e * sqrt(2.0 * pi()) * exp(0.5 * square(x));\r\n    x <- x - u / (1.0 + 0.5 * x * u);\r\n    return x;\r\n  }\r\n\r\n  vector invnorm(vector p, vector a, vector b, vector c, vector d, int N) {\r\n    real q; real r; real p_low;  real p_high; real e; real u;\r\n    vector[N] x;\r\n    p_low  <- 0.02425;\r\n    p_high <- 0.97575; # 1 - p_low    \r\n    for (i in 1:N) {\r\n      if( p[i] <= 0 )     x[i] <- negative_infinity();\r\n      else if (p[i] >= 1) x[i] <- positive_infinity();\r\n      else if ( (p_low <= p[i]) && (p[i] <= p_high) ) { # central region\r\n        q <- p[i] - 0.5;  \r\n        r <- square(q);\r\n        x[i] <- (((((a[1]*r+a[2])*r+a[3])*r+a[4])*r+a[5])*r+a[6])*q /\r\n                (((((b[1]*r+b[2])*r+b[3])*r+b[4])*r+b[5])*r+1.0);\r\n      }\r\n      else if (p[i] < p_low) { # lower region\r\n        q <- sqrt(-2.0*log(p[i]));\r\n        x[i] <- (((((c[1]*q+c[2])*q+c[3])*q+c[4])*q+c[5])*q+c[6]) /\r\n                ((((d[1]*q+d[2])*q+d[3])*q+d[4])*q+1.0);\r\n      }\r\n      else { # upper region\r\n        q <- sqrt(-2.0*log1m(p[i]));\r\n        x[i] <- -(((((c[1]*q+c[2])*q+c[3])*q+c[4])*q+c[5])*q+c[6]) /\r\n                 ((((d[1]*q+d[2])*q+d[3])*q+d[4])*q+1.0);\r\n      }\r\n      # refinement\r\n      e <- Phi(x[i]) - p[i];\r\n      u <- e * sqrt(2.0 * pi()) * exp(0.5 * square(x[i]));\r\n      x[i] <- x[i] - u / (1.0 + 0.5 * x[i] * u);\r\n    }\r\n    return x;\r\n  }\r\n\r\n  vector invnorm(real[] p, vector a, vector b, vector c, vector d, int N) {\r\n    return invnorm(to_vector(p), a, b, c, d, N);\r\n  }  \r\n}\r\n\r\ndata {\r\n  int<lower=1> N;\r\n}\r\ntransformed data {\r\n  vector[6] a; vector[5] b; vector[6] c; vector[4] d;\r\n  \r\n   a[1] <- -3.969683028665376e+01;\r\n   a[2] <-  2.209460984245205e+02;\r\n   a[3] <- -2.759285104469687e+02;\r\n   a[4] <-  1.383577518672690e+02;\r\n   a[5] <- -3.066479806614716e+01;\r\n   a[6] <-  2.506628277459239e+00;\r\n\r\n   b[1] <- -5.447609879822406e+01;\r\n   b[2] <-  1.615858368580409e+02;\r\n   b[3] <- -1.556989798598866e+02;\r\n   b[4] <-  6.680131188771972e+01;\r\n   b[5] <- -1.328068155288572e+01;\r\n\r\n   c[1] <- -7.784894002430293e-03;\r\n   c[2] <- -3.223964580411365e-01;\r\n   c[3] <- -2.400758277161838e+00;\r\n   c[4] <- -2.549732539343734e+00;\r\n   c[5] <-  4.374664141464968e+00;\r\n   c[6] <-  2.938163982698783e+00;\r\n\r\n   d[1] <-  7.784695709041462e-03;\r\n   d[2] <-  3.224671290700398e-01;\r\n   d[3] <-  2.445134137142996e+00;\r\n   d[4] <-  3.754408661907416e+00;\r\n}\r\nparameters {\r\n  real<lower=0,upper=1> p[N];\r\n}\r\nmodel {\r\n  // implicit i.i.d. uniform prior on p\r\n}\r\ngenerated quantities {\r\n  vector[N] x;\r\n  x <- invnorm(p, a, b, c, d, N); # should be ~ standard normal\r\n}\r\n\r\n```"},"comment":{"url":"https://api.github.com/repos/stan-dev/stan/issues/comments/77749383","html_url":"https://github.com/stan-dev/stan/issues/1157#issuecomment-77749383","issue_url":"https://api.github.com/repos/stan-dev/stan/issues/1157","id":77749383,"user":{"login":"hoanguc3m","id":7287528,"avatar_url":"https://avatars.githubusercontent.com/u/7287528?v=3","gravatar_id":"","url":"https://api.github.com/users/hoanguc3m","html_url":"https://github.com/hoanguc3m","followers_url":"https://api.github.com/users/hoanguc3m/followers","following_url":"https://api.github.com/users/hoanguc3m/following{/other_user}","gists_url":"https://api.github.com/users/hoanguc3m/gists{/gist_id}","starred_url":"https://api.github.com/users/hoanguc3m/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/hoanguc3m/subscriptions","organizations_url":"https://api.github.com/users/hoanguc3m/orgs","repos_url":"https://api.github.com/users/hoanguc3m/repos","events_url":"https://api.github.com/users/hoanguc3m/events{/privacy}","received_events_url":"https://api.github.com/users/hoanguc3m/received_events","type":"User","site_admin":false},"created_at":"2015-03-08T13:40:06Z","updated_at":"2015-03-08T13:40:06Z","body":"Is that what I should send? Rstan works fine in my computer. Thanks for your comments and instructions.\r\nhttps://github.com/hoanguc3m/Test-stan\r\n\r\n"}}