I am using  to JIT compile some looped python functions as part of a larger application. Ideally, everything will run in numba's "no python" mode, such that the loop can be parallelised.  One of these functions involves use of the Beta distribution CDF (the ):  $$I_z(z| \alpha, \beta) = \frac{B(z | \alpha, \beta)}{B(\alpha, \beta)}$$  where  $$B(\alpha, \beta) =\frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha + \beta)}$$  and  $$B(z|\alpha, \beta) = \int_0^{z} u ^{\alpha -1} (1 -u)^{\beta-1} du$$  To run in nopython mode, I would like to implement the above using only functions from python's standard math library and numpy.   Implementing $B(\alpha, \beta)$ is no problem as I can use math.gamma from the python standard library.   The problem is how to evaluate the integral for $B(z|\alpha, \beta)$, which -- based on my admittedly limited knowledge of the problem -- must be performed numerically.  Is there a "standard" way that this is implemented inside most statistical libraries, eg. scipy?  appears to ask a similar question but remains unanswered.  Please note that $\alpha$ and $\beta$ could be varied quite substantially, so a lookup table with interpolation is not an option.  Solution for those searching for similar: I ended up following  linked answer from the comment below by whuber. Doing some further reading lead me to  on Github, along with .
96327
Is there a way to closely approximate the CDF of a student-t distribution at a point $x$ without involving the hypergeometric function? For example, by using a series expansion, or expressing the CDF in terms of simpler functions?  The t-distribution can   $$t = \frac{X\sqrt{n}}{Y}$$  where $X$ is normally distributed with mean 0 and variance $\sigma^2$; $\frac{Y^2}{\sigma^2}$ has a chi-squared distribution. Maybe this expression can help simplify the t CDF?  I am trying to write an algorithm to compute the student-t CDF for any real number, so using the hypergeometric function will be extremely difficult. I am also trying to avoid integrating the PDF to get the CDF.
101315
I am looking at using the lasso as a method for selecting features and fitting a predictive model with a binary target.  Below is some code I was playing with to try out the method with regularized logistic regression.  My question is I get a group of "significant" variables but am I able to rank order these to estimate relative importance of each? Can the coefficients be standardized for this purpose of rank by absolute value (I understand that they are shown on the original variable scale through the coef function)? If so, how to do so (using the standard deviation of x and y) .  SAMPLE CODE:      library(glmnet)      #data comes from  #http://archive.ics.uci.edu/ml/datasets/Breast+Cancer+Wisconsin+(Diagnostic)      datasetTest &lt;- read.csv('C:/Documents and Settings/E997608/Desktop/wdbc.data.txt',head=FALSE)   #appears to use the first level as the target success    datasetTest$V2&lt;-as.factor(ifelse(as.character(datasetTest$V2)=="M","0","1"))   #cross validation to find optimal lambda #using the lasso because alpha=1      cv.result&lt;-cv.glmnet(                      x=as.matrix(dataset[,3:ncol(datasetTest)]),               y=datasetTest[,2],                       family="binomial",                       nfolds=10,                       type.measure="deviance",                      alpha=1                     )  #values of lambda used      histogram(cv.result$lambda)  #plot of the error measure (here was deviance) #as a CI from each of the 10 folds #for each value of lambda (log actually)      plot(cv.result)   #the mean cross validation error (one for each of the #100 values of lambda      cv.result$cvm  #the value of lambda that minimzes the error measure #result: 0.001909601      cv.result$lambda.min     log(cv.result$lambda.min)  #the value of lambda that minimzes the error measure #within 1 SE of the minimum #result: 0.007024236      cv.result$lambda.1se  #the full sequence was fit in the object called cv.result$glmnet.fit #this is same as a call to it directly. #here are the coefficients from the min lambda      coef(cv.result$glmnet.fit,s=cv.result$lambda.1se)
108794
I have no idea of Beta functions. Based on some properties of beta functions like, $B(m,n) = \int_{0}^\infty \frac{x^{m-1}}{(1+x)^{m+n}}dx \;\; ;B(m,n) = B(n,m)$  I arrived at $2B(m,n) = I+ \int_{1}^\infty \frac{x^{m-1}+x^{n-1}}{(1+x)^{m+n}}$ where I is my required integral. How to express the 2nd integral in terms of Beta function or is there any other way of solving. Need help. Thank You.
58580
I am trying to produce a model for which I have a response variable which is a proportion between 0 and 1, this includes quite a few 0s and 1s but also many values in between. I am thinking about attempting a beta regression. The package I have found for R (betareg) only allows values  in between 0 and 1 but not including 0 or 1 them selves. I have read elsewhere that theoretically the beta distribution should be able to handle values of 0 or 1 but I do not know how to handle this in R.I have seen some people add 0.001 to the zeros and take 0.001 from the ones, but I am not sure this is a good idea?  Alternatively I could logit transform the response variable and use linear regression. In this case I have the same problem with the 0 and 1's which cannot be log transformed.
43115
Suppose that  $X$ and $Y$ are two i.i.d. uniform random variables on the interval $[0,1]$  Let $Z=X/Y$, I am finding the cdf of $Z$, i.e. $ \Pr(Z\leq z) $.  Now, I came up with two ways of doing this. One produces a correct answer consistent with the pdf here: , the other does not. Why is the second method wrong?  First Method  $\newcommand{\rd}{\mathrm{d}} \Pr(Z\leq z) = \Pr(X/Y\leq z) = \Pr(X\leq zY) = \int^{1}_{0}\int^{\min(1,zy)}_{0} \rd x \rd y = \int^{1}_{0}\min(1,zy)\ \rd y$ $  = \left\{      \begin{array}{lr}        \int^{1/z}_{0}zy\ \rd y + \int^{1}_{1/z} \rd y&amp; : z &gt; 1\\        \int^{1}_{0}zy\ \rd y &amp; : z \leq 1      \end{array}    \right. $ $  = \left\{      \begin{array}{lr}        1 - \frac{1}{2z} &amp; : z &gt; 1\\        \frac{z}{2} &amp; : z \leq 1      \end{array}    \right. $  This appears correct.  Second Method  $\Pr(X/Y\leq z) = \Pr(X \leq zY\ |\ zY \geq 1)\Pr(zY \geq 1) + \Pr(X \leq zY\ |\ zY &lt; 1)\Pr(zY &lt; 1) $ by total probability  $ = \Pr(X \leq zY\ |\ zY \geq 1)\Pr(Y \geq 1/z) + \Pr(X \leq zY\ |\ zY &lt; 1)\Pr(Y &lt; 1/z)$  Taking $z&gt;1$ yields $ (1)(1-\frac{1}{z}) + (\int^{1/z}_{0}\int^{zy}_{0} \rd x \rd y)(\frac{1}{z}) = 1-\frac{1}{z} + (\int^{1/z}_{0}zy\ \rd y)(\frac{1}{z}) = 1-\frac{1}{z} + \frac{1}{2z^{2}}$  This is already different. Why is this wrong?  Thanks!
32918
According to Wikipedia, the  is related to the gamma distribution by the following relation:  $$\lim_{n\to\infty}n B(k, n) = \Gamma(k, 1)$$  Can you point me to a derivation of this fact? Can it be generalized? For example, is there a similar relation that results in something other than a constant 1 for the Gamma second parameter? What if we have  $$\lim_{n\to\infty,m\to\infty,n=mb}n B(k, m) $$  That is, the two variables go to infinity while maintaining a constant ratio b.  The reason I'm asking is because I'm trying to figure out how to simplify a hieraerchical bayesian model involving the beta distribution.  (This is my first post; sorry for the math notation, the MathJaX syntax was too daunting, but I'll try to learn)
22150
I have the following equations: \begin{equation} {1\over f(x,y)} {\partial f(x,y) \over \partial x} \alpha(x,y) + {1\over f(x,y)} {\partial f(x,y) \over \partial y} \beta(x,y) = \gamma(x,y) \end{equation} Where $\alpha(x,y)$, $\beta(x,y)$ and $\gamma(x,y)$ are exactly known functions. What I need to find is $f(x,y)$. Can anyone help me find a solution? Or are there any assumptions (other than $f(x,y) = g(x)h(y)$ that will allow me to solve this analytically? Is it a standard equation that I haven't come across? Many thanks guys!
58895
Been trying to this line to parse for a good 20mins, but all attempts have failed:  \beta = $-$0.84, t = $-$2.09, P = 0.04   Could anyone tell me why I'm getting a Missing $ inserted error here?
110641
I have a dynamic mixture distribution fitted to my risk data (i.e., I have all parameters) of Weibull and Generalized Pareto, with a Cauchy CDF mixing function, that can be written as:   \begin{align} \newcommand{\mixture}{{\rm mixture}} \newcommand{\Weibull}{{\rm Weibull}} \newcommand{\Cauchy}{{\rm Cauchy}} &amp;\mixture(x): x\in\mathbb{R^{+}} \to \mixture(x) \in [0,1]  \\ &amp;\mixture(x)=\big\{1-\Cauchy_{CDF}(x)\big\}\times\Weibull(x)+\Cauchy_{CDF}(x)\times GPD(x) \end{align}  in R, this gives:  mixture = function(x){   ((1 - pcauchy(x, location=535, scale=4.21e-04))*dweibull(x, shape=1.22, scale=62.31) +         pcauchy(x, location=535, scale=4.21e-04)*dgpd(x, xi=0.23, mu=0, beta=92.25))[1] }   I want to know if the following is the correct way of writing the quantile function of my mixture (where the $q$ subscripts stand for the quantile functions):  \begin{align} &amp;\mixture_{q}(y): y\in [0,1] \to \mixture_{q}(y) \in\mathbb{R^{+}}  \\ &amp;\mixture_{q}(y)=\big\{1-\Cauchy_{CDF}(\Weibull_{q}(y))\big\}\times \Weibull_{q}(y)+  \\ &amp;\hspace{32mm} \Cauchy_{CDF}(GPD_{q}(y))\times GPD_{q}(y) \end{align}  in R, this translates to:  mixture.quantile = function(y){   ((1 - pcauchy(qweibull(y, shape=1.22, scale=62.31), location=535, scale=4.21e-04)) *          qweibull(y, shape=1.22, scale=62.31) +          pcauchy(qgpd(y, xi=0.23, mu=0, beta=92.25), location=535, scale=4.21e-04) *          qgpd(y,xi=0.23,mu=0,beta=92.25))[1] }   The only thing that changes is that because we are going from $[0,1]$ to $\mathbb{R^{+}}$, but $\Cauchy_{CDF}$ goes from $\mathbb{R^{+}}$ to $[0,1]$, the weighing has to be performed based on the values of the quantile functions and cannot be done based on the values of $x$ anymore (hence the  $\Weibull_{q}(y)$ inside the $\Cauchy_{CDF}(.)$ for instance)... Is that correct??  With this quantile function I can get a Q-Q plot, or directly compare the quantile function of the mixture model with that of my data:  # load packages libs = c("repmis","fExtremes","evmix","evir") lapply(libs, library, character.only=TRUE)  # load data data = source_data("https://www.dropbox.com/s/r7i0ctl1czy481d/test.csv?dl=0")[,1]  # get quantile function values     mixture.quant = numeric(length=length(data)) vector.quant  = ppoints(length(data), a=0) for (i in 1:length(data)){   y                = vector.quant[i]   mixture.quant[i] = mixture.quantile(y) }  # Q-Q plot plot(mixture.quant, sort(data), xlab="theoretical quantile", ylab="sample quantile") abline(sort(data), sort(data), col="red")  # empirical CDF with mixture quantile function overlaid plot(mixture.quant, ppoints(length(data),a=0), col="red", type="l") lines(quantile(data, ppoints(length(data),a=0)), ppoints(length(data),a=0), type="l")       PS: see these other related threads for background: , .
114576
Compute: $$\lim_{x\to 0}\frac{\sqrt[m]{\cos(\alpha  x)}-\sqrt[n]{\cos{(\beta x)}}}{x^2},\;m,n\in\mathbb N,\alpha,\beta\in\mathbb R$$   My attempt:  $$L=\lim_{x\to 0}\frac{\sqrt[m]{\cos(\alpha x)}-\sqrt[n]{\cos{(\beta x)}}}{x^2}=\lim_{x\to 0}\frac{\sqrt[m]{\cos(\alpha x)}-1+1-\sqrt[n]{\cos(\beta x)}}{x^2}\\\displaystyle=\lim_{x\to 0}\frac{\left(1+(\cos(\alpha x)-1)\right)^{\frac{1}{m}}-1-\left(\left(1+(\cos(\beta x)-1)\right)^{\frac{1}{n}}-1\right)}{x^2}$$  Transformed expression: $$-\frac{(1+(\cos(\alpha x)-1))^{\frac{1}{m}}-1}{\cos(\alpha x)-1}\cdot\frac{1-\cos(\alpha x)}{\alpha^2x^2}\alpha^2+\frac{(1+(\cos(\beta x)-1))^{\frac{1}{n}}-1}{\cos(\beta x)-1}\cdot\frac{1-\cos(\beta x)}{\beta^2x^2}\beta^2$$  I applied the standard limits from the table: $$\lim_{x\to 0}\frac{(1+x)^a-1}{x}=a\;\&amp;\;\lim_{x\to 0}\frac{1-\cos x}{x^2}=\frac{1}{2}$$ $$L=-\frac{1}{m}\cdot\frac{1}{2}\alpha^2+\frac{1}{n}\cdot\frac{1}{2}\beta^2=\frac{m\beta^2-n\alpha^2}{2mn}$$ Is this correct?
108874
According to the , we have the scenario shown below, but how can ABC generate simulation datasets from samples of $\theta$ without knowing or evaluating the likelihood function?  For example, how would this method work for a normal $N(\mu, \sigma)$ distribution? Assuming that the parameter to estimate is $\mu$ and that $\sigma$ is known, how would it generate the simulations without evaluating $N(\mu, \sigma)$?
123715
I need to plot the Bessel functions of first and second kind (J and Y), and the modified Bessel functions of first and second kind (I and K), with integer order, from order 0 to order 5, using pgfplots. Several solutions have been already presented in other questions, but:  Gnuplot makes available for plots only the Bessel functions of the first and second kind and only of order 0 and 1, so  is not suitable. I would not like to use nor pstricks, neither LuaTeX, so also  and  are not suitable.  is quite similar, but has no answer and no examples in the comments.  How to accomplish this with pgfplots? Also an external solution with numpy, scipy and Matplotlib as in  would be ok (however, the code seems to work only for the J function).    I tried using the code in the last linked , inserting as an \addplot the file example-04.txt. It works for the function J, but it prints unacceptable values (e+09) for function Y, which is very large and negative near 0. The code which generates example-04.txt is Python and here it would be almost certainly OT.
15933
I've been thinking about this problem.  The usual logistic function for modeling binary data is: $$ \log\left(\frac{p}{1-p}\right)=\beta_0+\beta_1X_1+\beta_2X_2+\ldots $$ However is the logit function, which is an S-shaped curve, always the best for modeling the data?  Maybe you have reason to believe your data does not follow the normal S-shaped curve but a different type of curve with domain $(0,1)$.  Is there any research into this?  Maybe you can model it as a probit function or something similar, but what if it is something else entirely?  Could this lead to better estimation of the effects?  Just a thought I had, and I wonder if there is any research into this.
120164
Find the ordered pair $(\alpha,\beta)$ with non-infinite $\beta \ne 0$ such that $$\lim_{n\to\infty} \frac{\sqrt[n^2]{1!2!\cdots n!}}{n^\alpha} = \beta$$  My approach:   $$\ln (1!2!\cdots n!) = (n)\ln 1 + (n-1)\ln 2 + \cdots + (2)\ln (n-1) + \ln(n) \\ \begin{align} = n\ln\left(\frac{1}{n}\right) + (n-1)\ln\left(\frac{2}{n}\right) + \cdots + \ln\left(\frac{n}{n}\right) + \frac{(n)(n+1)}{2} \ln (n)\end{align}$$  Then $$\ln(\sqrt[n^2]{1!2!\cdots n!}) = \frac{1}{n^2} \ln(1!2!\cdots n!) = \frac{n+1}{n} \cdot \ln n + \frac{1}{n} \left[\sum_{m=1}^n \left(\frac{n+1-m}{n} \cdot \ln \frac{m}{n}\right)\right]$$  And that's about as far as I got. Any ideas about proceeding with this method or perhaps even with a different method?  Thanks  A
62583
I have a dataset as follows:  ni-xi 76 48 38 19 20 14 26 29 39 45 38 36 32 34 26 21 23 12 24 18 14 61 35 21 32 20 89 30 34 35  xi    0  0 16 26  1  8 23 19 24 16 16 19  0  0  0  6  8  3  5  9  5  0 34 17 13  0  1 11 13 22   I want to model this dataset as a beta-binomial distribution: $x_i \sim \mathcal{Binomial}(p, n_i)$ where $p \sim \mathcal{Beta}(α, β)$.  I use VGAM package to do the fitting.   fit=vglm(cbind(ni-xi, xi) ~ 1, betabinomial.ab, trace = TRUE) coef(fit)  VGLM    linear loop  1 :  loglikelihood = -95.9867 VGLM    linear loop  2 :  loglikelihood = -95.6633 VGLM    linear loop  3 :  loglikelihood = -95.6483 VGLM    linear loop  4 :  loglikelihood = -95.6477 VGLM    linear loop  5 :  loglikelihood = -95.6477 VGLM    linear loop  6 :  loglikelihood = -95.6477   However, I get the results as follows:  (Intercept):1 (Intercept):2      0.9006006    -0.3794863   Can anyone tell me why I get the negative parameters? How can i solve this problem.
15769
The  is a discrete probability distribution over a set $\{1,\dots, N\}$ with the probability mass function  $$ p(1)=\frac{1}{N},\qquad p(k)=\frac{1}{k(k-1)}\quad\text{for }k\in\{2,\dots, N\} $$  I'd like to use it as part of an implementation of an , ideally in Python where a uniform random number generator is available.
22014
I was doing my homework and the following problem came up! We have the CDF of an exponential distribution that is shifted $L$ units where $L&gt;0$ and $x&gt;=L$. The CDF is:  $$1-e^{-\lambda(x-L)}$$  The question says that we should assume that the following data are lifetimes of electric motors, in hours, which are:  $$\begin{align*} 153.52,103.23,31.75,28.91,37.91,7.11,99.21,31.77,11.01,217.40 \end{align*}$$  Please note that the $mean$ of these numbers is: $72.182$  Now the question has two parts which I will go through one by one:  Part1: Evaluate the log likelihood for the data when $\lambda=0.02$ and $L=3.555$. Now the way I approached the problem was to take the derivative of the CDF with respect to $\lambda$ to get the PDF which is:  $$(x-L)e^{-\lambda(x-L)}$$  Then since we have $n$ observations where $n=10$, we have the following joint pdf, due to independence:  $$(x_i-L)^ne^{-\lambda(x_i-L)n}$$ which can be rewritten as the following log likelihood:  $$n\ln(x_i-L)-\lambda\sum_{i=1}^n(x_i-L)$$ Is this the correct approach? Because it would take quite a while and be pretty cumbersome to evaluate $n\ln(x_i-L)$ for every observation?  Part2: The question also asks for the ML Estimate of $L$. So assuming the log likelihood is correct, we can take the derivative with respect to $L$ and get:  $\frac{n}{x_i-L}+\lambda=0$ and solve for $L$? Is this correct? Because I am not quite sure on how I should proceed?  Thanks so much for your help! I greatly appreciate it :)