// Mata functions NJC 15 December 2006 mata : // vector = f(y | alpha, beta) real invgammaden(real scalar alpha, real scalar beta, real colvector y) { if (alpha <= 0) { displayas("error") printf("alpha must be positive\n") exit(error(411)) } if (beta <= 0) { displayas("error") printf("beta must be positive\n") exit(error(411)) } if (sum(y :<= 0)) { displayas("error") printf("values must be positive\n") exit(error(411)) } return((y:^(-(alpha + 1)) :* exp(- beta :/ y)) / (beta^(-alpha) * gamma(alpha))) } // variable = f(y | alpha, beta) void invgammadenvar(real scalar alpha, real scalar beta, string scalar varname, string scalar tousename, string scalar densityname) { real colvector y, density if (alpha <= 0) { displayas("error") printf("alpha must be positive\n") exit(error(411)) } if (beta <= 0) { displayas("error") printf("beta must be positive\n") exit(error(411)) } y = st_data(., varname, tousename) if (sum(y :<= 0)) { displayas("error") printf("values must be positive\n") exit(error(411)) } density = (y:^(-(alpha + 1)) :* exp(- beta :/ y)) / (beta^(-alpha) * gamma(alpha)) st_addvar("double", densityname) st_store(., densityname, tousename, density) } // vector = F(y | alpha, beta) real invgamma(real scalar alpha, real scalar beta, real colvector y) { real colvector dist if (alpha <= 0) { displayas("error") printf("alpha must be positive\n") exit(error(411)) } if (beta <= 0) { displayas("error") printf("beta must be positive\n") exit(error(411)) } if (sum(y :<= 0)) { displayas("error") printf("values must be positive\n") exit(error(411)) } dist = 1 :- gammap(alpha, beta :/ y) return(dist) } // variable = F(y | alpha, beta) void invgammavar(real scalar alpha, real scalar beta, string scalar varname, string scalar tousename, string scalar distname) { real colvector y, dist if (alpha <= 0) { displayas("error") printf("alpha must be positive\n") exit(error(411)) } if (beta <= 0) { displayas("error") printf("beta must be positive\n") exit(error(411)) } y = st_data(., varname, tousename) if (sum(y :<= 0)) { displayas("error") printf("values must be positive\n") exit(error(411)) } dist = 1 :- gammap(alpha, beta :/ y) st_addvar("double", distname) st_store(., distname, tousename, dist) } // vector = quantile(P | alpha, beta) real invinvgamma(real scalar alpha, real scalar beta, real colvector p) { real colvector y if (alpha <= 0) { displayas("error") printf("alpha must be positive\n") exit(error(411)) } if (beta <= 0) { displayas("error") printf("beta must be positive\n") exit(error(411)) } if (sum(p :<= 0 :| p :>= 1)) { displayas("error") printf("values must be in (0,1)\n") exit(error(498)) } y = beta :/ invgammap(alpha, 1 :- p) return(y) } // variable = quantile(P | alpha, beta) void invinvgammavar(real scalar alpha, real scalar beta, string scalar varname, string scalar tousename, string scalar yname) { real colvector p, y if (alpha <= 0) { displayas("error") printf("alpha must be positive\n") exit(error(411)) } if (beta <= 0) { displayas("error") printf("beta must be positive\n") exit(error(411)) } p = st_data(., varname, tousename) if (sum(p :<= 0 :| p :>= 1)) { displayas("error") printf("values must be in (0,1)\n") exit(error(498)) } y = beta :/ invgammap(alpha, 1 :- p) st_addvar("double", yname) st_store(., yname, tousename, y) } end