# student t pdf

Here's a code for student t pdf.

var("s n t")
assume(n, 'integer')
assume(n>0)
tdist=integral(sqrt(s/n)*1/sqrt(2*pi)*exp(-s*t^2/(2*n))*exp(-s/2)*s^(n/2-1)/(gamma(n/2)*2^(n/2)),s,0,infinity)
tdist


And it says, TypeError: Computation failed since Maxima requested additional constraints (try the command 'assume(n+1)/2>0)' before integral or limit evaluation, for example): Is (n+1)/2 an integer?

The funny thing is when I add n=5 or n=3, there's no error!

var("s n t")
assume(n, 'integer')
n=3
tdist=integral(sqrt(s/n)*1/sqrt(2*pi)*exp(-s*t^2/(2*n))*exp(-s/2)*s^(n/2-1)/(gamma(n/2)*2^(n/2)),s,0,infinity)
tdist


EDIT : Actually, I tried this on MATHLAB to fail.

   My last hope would be Mathematica? wolframalpha.com says it's timed out!

edit retag close merge delete

I think that once you say n=3, it's sort of substituted in. If you input sqrt(s/n)*1/sqrt(2*pi)*exp(-s*t^2/(2*n))*exp(-s/2)*s^(n/2-1)/(gamma(n/2)*2^(n/2)) you get 1/2*sqrt(1/3)*s*e^(-1/6*s*t^2 - 1/2*s)/pi, no n for Maxima to worry about.

( 2012-06-07 16:06:45 +0100 )edit

The problem is I want the integration for general n(positive integer), not for some specific n=3 or 5!

( 2012-06-08 03:06:19 +0100 )edit

I understand, I was just pointing out why there was no error.

( 2012-06-08 10:32:46 +0100 )edit

Sort by ยป oldest newest most voted

I don't think we're going to be able to fix this very easily at all. The problem is Maxima's framework, which there has been much discussion about and which is not going to change.

(%i5) assume(n>0);
(%o5)                               [n > 0]
(%i6) display2d:false;

(%o6) false
(%i7) integrate(sqrt(s/n)*1/sqrt(2*pi)*exp(-s*t^2/(2*n))*exp(-s/2)*s^(n/2-1)/(gamma(n/2)*2^(n/2)),s,0,inf);

Is (n+1)/2 an integer?

no
;
(%o7) gamma(n/2+1/2)*2^(-n/2-1/2)*(t^2/(2*n)+1/2)^(-n/2-1/2)
/(sqrt(n)*gamma(n/2)*sqrt(pi))
(%i8) integrate(sqrt(s/n)*1/sqrt(2*pi)*exp(-s*t^2/(2*n))*exp(-s/2)*s^(n/2-1)/(gamma(n/2)*2^(n/2)),s,0,inf);

Is (n+1)/2 an integer?

yes;
Is n an integer?

yes;
(%o8) gamma(n/2+1/2)*2^(-n/2-1/2)*(t^2/(2*n)+1/2)^(-n/2-1/2)
/(sqrt(n)*gamma(n/2)*sqrt(pi))


Don't even ask about why it asks if n is an integer when it knows (n+1)/2 is an integer; Maxima itself acknowledges that its assumption framework is weak, and of course every open source project has limited resources to fix big architectural issues like that. I don't think it would be trivial for them to fix that.

As for Sage, assuming n is an integer will help a little, but only so much.

(%i9) declare(n,integer);

(%o9) done
(%i10) integrate(sqrt(s/n)*1/sqrt(2*pi)*exp(-s*t^2/(2*n))*exp(-s/2)*s^(n/2-1)/(gamma(n/2)*2^(n/2)),s,0,inf);

Is (n+1)/2 an integer?

yes;
(%o10) gamma(n/2+1/2)*2^(-n/2-1/2)*(t^2/(2*n)+1/2)^(-n/2-1/2)
/(sqrt(n)*gamma(n/2)*sqrt(pi))

(%i11) declare((n+1)/2,integer);

declare: improper argument: (1+n)/2
-- an error. To debug this try: debugmode(true);


So we wouldn't even be able to try that, even if we wanted to; Maxima only lets "atoms" be declared to be something (i.e., variables), though we can assume positive or negative about them. And even though the answers are the same, many have pointed out that going through a decision tree of this kind and trying to check that can be arbitrarily computationally expensive, if it's even possible.

I'm sorry, but this is something at the core of how Maxima does business. On the plus side, the answer works nicely enough that one gets from Maxima, and you could use it in Sage.

more

What's about substituting n to make it an odd integer?

var("s k n t")
assume(k, 'integer')
assume(k>0)
term= sqrt(s/n)*1/sqrt(2*pi)*exp(-s*t^2/(2*n))*exp(-s/2)*s^(n/2-1)/(gamma(n/2)*2^(n/2))
tdist=integral(term.substitute(n=2*k+1),s,0,infinity)
tdist

more

Yes, it works! But what about when n is an even integer? It also has solution when n is an even number, which is of course p.d.f of t-distribution!

( 2012-06-07 08:08:11 +0100 )edit

Does wrapping the expression with a python function fit your needs?

var("s n t")
term= sqrt(s/n)*1/sqrt(2*pi)*exp(-s*t^2/(2*n))*exp(-s/2)*s^(n/2-1)/(gamma(n/2)*2^(n/2))

def tdist(k):
return integral(term.substitute(n=k),s,0,infinity)

tdist(4)

more

Thank you for the answer. But the full solution should be f(t,n). we could consider f(t,n) from f(t,1), f(t,2), f(t,3)... but it's not the exact solution!!!

( 2012-06-08 03:01:49 +0100 )edit
var("s t")
def f(t,n):
return integral(sqrt(s/n)*1/sqrt(2*pi)*exp(-s*t^2/(2*n))*exp(-s/2)*s^(n/2-1)/(gamma(n/2)*2^(n/2)),s,0,infinity)

print [f(t,k) for k in range(1,10)]

more

This is a modification of the solution by ndomes

( 2012-06-08 04:34:02 +0100 )edit

On the off chance you just need the Student t distribution:

See this Sage documentation and search for t-distribution.

sage: nu = 1
sage: T = RealDistribution('t', nu)
sage: T.get_random_element() # random
-0.994514581164
sage: T.distribution_function(0)
0.318309886184
sage: T.cum_distribution_function(1)
0.75
sage: T.cum_distribution_function_inv(.5)
0.0


The word "student's" was removed in the change associated with Trac 9080 which added F-distribution. It will be a lot more work to have the distributions well integrated with the rest of Sage, but we could put that word back in - I suppose it was my fault, since I reviewed the ticket.

Also, see this sagenb worksheet for someone plotting the actual pdf.

more