Ask Your Question
0

student t pdf

asked 2012-06-07 06:44:14 +0100

anonymous user

Anonymous

updated 2022-09-01 08:17:52 +0100

FrédéricC gravatar image

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

Could any one help me about this?

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

   My last hope would be Mathematica? wolframalpha.com says it's timed out!
edit retag flag offensive close merge delete

Comments

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.

kcrisman gravatar imagekcrisman ( 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!

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

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

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

5 Answers

Sort by » oldest newest most voted
1

answered 2012-06-08 10:54:06 +0100

kcrisman gravatar image

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.

edit flag offensive delete link more
1

answered 2012-06-07 07:23:30 +0100

ndomes gravatar image

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
edit flag offensive delete link more

Comments

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!

QuantLabId gravatar imageQuantLabId ( 2012-06-07 08:08:11 +0100 )edit
0

answered 2012-06-08 04:32:02 +0100

achrzesz gravatar image
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)]
edit flag offensive delete link more

Comments

This is a modification of the solution by ndomes

achrzesz gravatar imageachrzesz ( 2012-06-08 04:34:02 +0100 )edit
0

answered 2012-06-07 13:55:34 +0100

ndomes gravatar image

updated 2012-06-08 03:39:19 +0100

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)
edit flag offensive delete link more

Comments

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!!!

QuantLabId gravatar imageQuantLabId ( 2012-06-08 03:01:49 +0100 )edit
0

answered 2012-06-08 11:00:01 +0100

kcrisman gravatar image

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.

edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2012-06-07 06:44:14 +0100

Seen: 800 times

Last updated: Jun 08 '12